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Abstract 

Relative  orbit  elements  (ROEs)  based  on  a  circular  chief  satellite  orbit  are  er¬ 
roneous  when  applied  to  a  perturbed,  non-circular  reference  orbit.  In  those  situa¬ 
tions,  the  ROEs  will  encounter  geometric  instability  and  drift.  To  counter  this,  a 
set  of  time-variant  ROEs  have  been  derived  to  describe  the  relative  orbit  for  both 
the  unperturbed,  elliptical  chief  and  the  perturbed,  circular  chief.  A  highly  coupled 
relationship  is  found  that  describes  the  relative  trajectory  to  higher  accuracy  when 
compared  to  numeric  integration.  To  show  the  applicability  of  the  ROEs  to  formation 
design,  methods  to  initialize  a  stationary  relative  orbit  are  detailed  and  an  algorithm 
for  ROE-based  guidance  and  navigation  is  proposed.  The  results  provide  a  method 
to  predict  the  relative  motion,  while  examining  time-varying  parameters  of  the  mo¬ 
tion.  Eccentricity  effects  are  shown  to  induce  severe  time-variance  to  the  system  and 
introduce  a  level  of  mathematical  abstraction.  Perturbing  J2  effects  are  shown  to 
introduce  periodic  effects  and  compound  the  secular  variations  to  the  circular  ROEs. 
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Osculating  Relative  Orbit  Elements  Resulting  erom 
Chiee  Eccentricity  and  J2  Perturbing  Eorces 

I.  Introduction 

Relative  satellite  motion  has  seen  vast  usage  and  progress  in  recent  years.  Nu¬ 
merous  well-defined  solution  sets  mapping  the  motion  of  a  reference  satellite,  the 
chief,  with  a  target  satellite,  the  deputy,  have  been  developed.  Dynamic  models  of 
satellite  motion  strongly  influence  close  proximity  operations  including  rendezvous 
and  docking.  A  past  example  of  proximity  operations  exists  in  the  Apollo  era,  during 
which  dynamic  algorithms  were  utilized  for  modular  rendezvous  [1].  However,  these 
operations  were  conducted  under  the  assistance  of  astronauts  in  the  loop,  not  allow¬ 
ing  for  full  autonomy.  The  idea  of  autonomous  close  proximity  operations  correlates 
with  and  motivates  the  attempt  to  further  generalize  the  dynamics  of  relative  satellite 
motion. 

As  a  vital  role  in  the  Apollo  program,  these  special  operations  were  studied 
primarily  for  docking  procedures.  Current  operations  still  maintain  focus  on  the 
same  docking  and  rendezvous  problem.  In  order  to  actively  control  this  maneuver,  the 
relative  separation  between  the  two  bodies  is  gradually  forced  to  zero.  For  example, 
the  resupply,  modification,  and  construction  of  the  International  Space  Station  is 
effected  by  the  use  of  autonomous  unmanned  orbit  transfer  vehicles,  relying  heavily 
on  control  of  the  rendezvous  problem  [2]. 

Formation  flight  is  another  defining  example  of  the  application  of  dynamic  mod¬ 
els  for  relative  motion  between  two  or  more  satellites.  By  utilizing  several  simple  and 
smaller  satellites,  the  size  and  complexity  of  single  large  spacecraft  missions  can  be 
reduced.  For  example,  a  reconnaissance  satellite  scaled  to  support  a  large  aperture 
could  be  reduced  to  a  formation  of  multiple  smaller  bodies  to  support  the  same  mis¬ 
sion  [3].  Failsafe  options  are  also  inherent  in  formation  flight,  in  that  a  catastrophic 
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failure  in  one  subsystem  does  not  necessarily  imply  mission  failure;  this  is  known  as 
’’gradual  degradation.”  Dynamic  models  of  relative  motion  can  be  applied  to  for¬ 
mations  by  selecting  a  reference  satellite  and  propagating  the  relative  equations  of 
motion  simultaneously  for  each  deputy,  setting  a  design  point  via  initial  conditions. 

Any  model  used  in  real-time  application  will  tend  towards  a  trade  between 
complexity  and  accuracy.  The  representation  of  every  possible  force  on  a  body  is 
impossible,  and  also  unnecessary  in  practice.  The  model  can  be  looked  at  in  terms  of 
kinematic  and  kinetic  terms.  The  kinematic  terms  will  include  the  differential  eqna- 
tions  mapping  the  motion,  while  the  kinetic  terms  can  be  thonght  of  as  forcing  terms 
to  the  motion.  From  a  deterministic  standpoint,  the  kinematic  part  of  the  model  can 
be  found,  while  the  forcing  terms  are  set  and  assnmed  for  the  application.  More  than 
likely,  the  forcing  terms  will  be  analytically  expressed  throngh  kinematic  variables, 
implying  that  the  inclusion  of  various  forces  will  directly  impact  the  coupling  and 
complexity  of  the  differential  eqnations  of  motion  for  the  model.  The  assnmptions 
made  in  each  model  will  have  a  substantial  mission  impact.  Software  constrncted 
in  accordance  with  analytical  Newtonian  dynamics  can  often  be  implemented  in  an 
autonomous  loop.  In  this  regard,  there  is  trade  off  between  model  simplicity  and 
dynamical  error  in  the  software.  In  addition,  the  simplicity  of  the  model  drives  the 
computing  power  needed  to  numerically  solve  the  differential  equations.  This  would 
not  be  a  problem  if  a  generalized  analytical  solntion  of  relative  satellite  motion  conld 
be  derived;  however,  the  stochastic  pertnrbations  of  spaceflight  allow  only  for  approx¬ 
imations  and  nnmerical  solutions  based  on  a  priori  operations. 

Inherent  in  these  models  are  assnmptions  constraining  the  orbit  of  the  chief,  and, 
in  addition,  the  proximity  between  the  two  satellites.  The  most  common  assumptions 
are  those  reqnired  for  the  linearization  of  the  model  derived  by  Clohessy  and  Wilt¬ 
shire  [4]  (HCW):  a  circular  chief,  unperturbed  two-body  Keplerian  dynamics,  and 
close  proximity.  The  HCW  model  often  manifests  itself  with  respect  to  rendezvous, 
docking  manenvers,  and  the  dynamics  of  satellite  clusters  [5],  but  the  linearization 
carries  with  it  limitations  in  real  world  applications. 
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Of  high  importance  to  the  present  study  is  one  particular  realization  of  the  HCW 
solution,  which  is  the  parameterization  of  the  resulting  trajectory  by  six  relative  orbit 
elements  (ROEs).  The  HCW  model  will  be  discussed  in  great  detail  later  in  this 
study,  along  with  the  derivation  of  the  ROEs  and  their  utility. 

The  linear  time-invariant  ordinary  differential  equations  of  the  HCW  formula¬ 
tion  allow  for  an  elegant  closed  form  solution.  A  fallout  of  the  HCW  model  shows 
that  the  relative  orbit  between  the  deputy  and  the  chief  will  be  an  ellipse  with  a  semi¬ 
major  axis  that  is  twice  the  semi-minor  axis  (a  2x1  ellipse)  when  projected  into  the 
orbital  plane  of  the  chief.  The  previously  mentioned  ROEs  provide  an  instantaneous 
snapshot  of  the  relative  orbit  that  describes  the  size,  shape,  and  orientation  of  the 
projected  relative  ellipse.  The  derivation  of  the  ROEs  is  performed  under  the  iden¬ 
tical  assumptions  as  the  HCW  model.  The  inclusion  of  orbital  perturbations  to  the 
HCW  model  represents  an  insertion  of  non-linearity  and  time-variance  to  the  ROEs. 
In  that  the  ROEs  are  linear  and  non-linear  combinations  of  the  relative  state  vector 
and  a  stand-alone  parameterization,  the  question  stands  as  to  what  impacts  are  made 
on  the  ROEs  when  operating  with  a  perturbed  or  non-circular  chief.  These  resulting 
relative  orbit  elements  are,  in  the  context  of  this  study,  defined  as  osculating  relative 
orbit  elements.  To  increase  the  utility  and  applicability  of  the  relative  orbit  element 
parameterization,  the  osculating  relative  orbit  response  to  both  the  relaxation  of  the 
unperturbed  motion  assumption  and  the  dynamics  of  the  non-circular  chief  is  the  aim 
of  the  current  study.  The  following  section  provides  detail  on  the  individual  aspects 
of  the  study. 

1.1  Problem  Statement 

The  purpose  of  this  study  is  to  develop  numerical  and  analytical  results  for 
perturbations  in  the  relative  orbit  elements  initiated  by  chief  eccentricity  and  the  J2 
zonal  harmonic.  The  study  will  be  divided  into  two  primary  phases,  investigating  the 
eccentricity  and  J2  effects  separately,  followed  by  proposed  applications. 
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1.1.1  Phase  One:  Osculating  ROEs  Perturbed  by  Chief  Eccentricity.  The 
first  phase  will  be  a  look  at  how  the  relative  orbit  elements  vary  as  a  result  of  a 
non-circular  chief.  Numerically  integrated  two-body  Keplerian  motion  (2BP)  will  be 
assumed  as  truth,  with  the  2BP  state  vector  numerically  substituted  as  the  arguments 
for  the  ROE  analytic  functions.  Stationary  orbits  will  be  primarily  investigated. 
Using  the  Yamanaka-Ankersen  (YA  [6])  state  transition  matrix  for  relative  motion, 
full  expressions  for  the  ROEs  will  be  derived  as  a  function  of  initial  conditions.  The 
objective  of  this  exercise  is  to  determine  the  validity  and  applicability  of  ROEs  for 
geometric  visualization  of  the  unperturbed  relative  orbit. 

1.1.2  Phase  Two:  Osculating  ROEs  Perturbed  by  the  .J2  Harmonic.  Ignoring 
atmospheric  drag,  the  spherical  central  body  assumption  in  the  Clohessy-Wiltshire 
equations  serves  as  a  primary  error  source.  The  second  phase,  following  in  a  similar 
fashion  as  Phase  One,  will  examine  the  behavior  of  the  ROEs  for  a  circular  chief 
perturbed  by  the  J2  harmonic.  Having  been  well  documented  in  the  literature  as  a 
well-behaved  model  for  the  J2  perturbed  chief,  the  Schweighart-Sedwick  [7]  model  is 
fully  utilized  to  derive  analytical  expressions  for  the  ROEs  as  a  function  of  initial 
conditions.  The  Gim-Alfriend  [8]  model  will  be  used  as  a  purely  analytical  study  to 
understand  the  physical  effects.  Being  that  the  ROEs  are  a  HOW  realization,  baseline 
comparisons  are  made  to  the  unperturbed  HOW  model. 

1.1.3  Applications  of  Osculating  ROEs. 

1.1. 3.1  Stationary  Orbit  Initialization.  To  provide  a  means  of  for¬ 
mation  installation  utilizing  ROEs,  sets  of  initial  conditions  to  allow  for  stationary 
relative  orbits  is  derived  using  period  matching  and  methods  of  J2  invariance. 

1.1. 3. 2  Guidance  and  Maneuver.  As  an  application  of  the  ROEs, 
a  guidance  algorithm  is  proposed  to  modify  the  orbit  using  osculating  ROEs.  The 
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purpose  of  this  study  is  to  introduce  the  concept  and  applicability  of  controls  to 
osculating  ROEs. 
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II.  Background 

The  relative  motion  between  two  bodies  is  not  a  new  concept.  More  than  likely,  the 
mapping  of  the  motion  of  one  object  compared  to  that  of  another  is  a  daily  occurance. 
Walking  side  by  side,  driving  on  a  highway,  running  along  a  track;  the  examples  are 
endless.  Luckily,  rather  than  free  will  driving  the  motion  of  one  body,  the  differential 
equations  of  orbital  mechanics  determine  the  relative  trajectories.  While  the  assump¬ 
tion  of  Keplerian  motion  greatly  simplihes  the  governing  laws,  various  perturbations 
can  increase  the  complexity  of  the  dynamical  models.  This  background  will  provide 
an  overview  of  the  past  and  current  state  of  the  art.  Section  2.1  will  provide  historical 
examples  of  relative  motion  applications.  Sections  2.2  to  2.4  will  provide  the  basic 
background  and  reference  frames  needed  to  understand  the  mathematics  in  the  clas¬ 
sical  ROE  expressions.  Finally,  Section  2.5  will  discuss  the  impact  of  the  non-circular 
chief  and  the  J2  perturbation  on  relative  satellite  motion. 

2.1  Relative  Satellite  Motion  Throughout  History 

2.1.1  Historical  Examples.  Between  June  of  1983  and  August  of  2005,  a 
total  of  57  shuttle  missions  utilized  one  or  more  forms  of  close  proximity  operations 
successfully  [9].  The  objectives  of  these  operations  vary  from  formation  maintenance 
to  hnal  docking  with  a  desired  chief.  Prior  to  these  maneuvers,  experiments  to  validate 
the  ability  of  a  human  eye  to  track  and  maintain  manual  control  of  a  docking  sequence 
were  performed  on  Mercury  missions.  Following  the  Mercury  Program,  Gemini  set 
forth  and  established  a  solid  foundation  for  the  future  of  human  rendezvous.  Good¬ 
man  states  that  the  most  significant  accomplishments  of  the  Gemini  Program  with 
respect  to  rendezvous  operations  included  high  and  low  orbit  coelliptic  rendezvous,  or¬ 
bital  night  and  day  docking,  optical  measurement  rendezvous,  conjunctive  countdown 
for  the  maneuver  from  both  chief  and  deputy  perspective,  and  multiple  rendezvous 
operations  while  staying  in  a  propellant  budget  [9] .  The  Apollo  Program  established 
rendezvous  operations  as  methodical  techniques,  using  several  missions  to  practice  lu¬ 
nar  landing.  As  a  side  note,  the  idea  of  lunar  rendezvous  was  not  a  new  thought  at  this 
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time;  one  of  the  first  accounts  of  relative  motion  model  development  is  Hill’s  [10]  lunar 
equations:  his  attempt  to  describe  relative  motion  in  the  Earth-Moon-Sun  system. 

Moving  the  focus  to  that  of  formation  flying,  the  size,  cost,  and  complexity 
of  many  classes  of  missions  can  be  greatly  reduced  by  utilizing  a  number  of  smaller 
satellites.  Using  multiple  satellites  allows  multi-tasking  within  the  formation.  In  the 
event  of  a  catostrophic  failure  in  one  of  the  deputy  satellites  within  the  formation, 
adjustments  can  be  made  to  compensate  for  the  loss.  [3].  Examples  of  formation 
flight  include  the  Orion  program  as  a  proof  of  concept  for  GPS  based  relative  nav¬ 
igation,  and  the  ESA  Cluster  mission  for  magnetospheric  studies.  Stationkeeping 
requirements  for  formation  control  tend  to  be  on  a  low  order  of  required  thrust.  This 
has  led  to  direct  development  of  electric  thrusters  utilizing  electromagnetic  as  well  as 
electrostatic  forces  [11].  A  formation  design  will  specify  initial  conditions  in  such  a 
way  to  satisfy  the  mission  requirements.  For  example,  a  circular  in-plane  formation 
will  be  bounded  by  specifying  the  in-track  and  cross-track  states  as  a  function  of  the 
radial  states.  While  specifying  these  values  limits  the  the  initial  degrees  of  freedom  (in 
this  case,  only  having  two  initial  conditons  to  choose),  the  motion  is  deterministically 
known  and  the  free  initial  conditions  can  be  used  to  set  the  period  and  phasing  of  the 
deputy  [11]. 

2.1.2  Model  Development.  Whether  having  been  applied  to  rendezvous 
operations  or  to  the  control  of  a  formation,  the  idea  of  relative  motion  models  has 
been  widely  used  throughout  the  history  of  the  space  program.  The  history  of  model 
development  in  this  held  is  an  interesting  study  in  various  attempts  at  hrst  order 
linearizations,  linearized  state  transition  matrices,  incorporating  various  orbital  per¬ 
turbations,  and  increasing  model  hdelity.  In  fact,  the  history  of  model  development 
can  almost  be  seen  as  a  survey  of  the  state  transition  matrices  derived  over  the  years. 
The  idea  of  propagation  of  an  initial  state  to  a  desired  time  through  determinstic 
dynamics  is  quite  popular  in  both  theory  and  practice.  The  popularity  is  of  course 
due  to  accuracy  without  time-consuming  and  computationally  expensive  numerical 
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integration.  The  initiation  of  a  state  transition  matrix  depends  on  the  linearization 
inherent  in  the  model.  Moreover,  the  independent  variable  used  in  the  propagation 
is  of  extreme  importance.  For  example,  if  time  is  used  then  certain  states  may  be 
poorly  represented  near  perigee  without  the  appropriate  time  step;  however,  if  true 
anomaly  is  used,  complicating  dehnite  integrals  will  show  themselves.  A  brief  survey 
on  the  properties  of  state  transition  matrices  can  be  found  in  Appendix  A.  Despite 
the  broad  dehnitions  of  these  models.  Carter  [12]  proposes  three  classihcations  of 
linearized  models 

1.  Inertial  or  rotating  reference  frames 

2.  Linearizing  the  central  force  held  to  introduce  a  gradient  to  the  equations  of 
motion,  or  linearizing  orbital  parameters 

3.  Nominal  reference  orbit 

Although  proving  more  abstract  to  view  initially,  the  rotating  frame  is  employed 
the  most  often  as  a  means  to  visualize  system  behavior  in  the  vicinity  of  the  chief. 
More  specihcally,  the  linearization  of  an  inverse  square  gravity  term  (p/r^)  yields 
the  most  well  known  relative  motion  model  in  the  held:  Clohessy  and  Wiltshire  [4], 
which  is  discussed  in  great  detail  in  Section  2.3.  Limitations  of  this  model  are  a 
result  of  the  linearizations  made  to  derive  the  closed  form  solution.  Constraints  such 
as  a  circular  chief,  proximity  assumptions,  and  unperturbed  motion  limit  the  long 
term  applications  of  the  model.  The  circular  chief  was  expanded  to  second-order 
truncation  by  Karlgaard  and  Lutze  using  the  method  of  multiple  scales  in  spherical 
coordinates  [13]  with  signihcant  increase  in  accuracy  over  HCW. 

The  circular  chief  assumption  is  highly  limiting,  most  paramount  in  that  a  pure 
circular  orbit  is  not  plausible.  For  missions  involving  vehicle  servicing  or  sample  re¬ 
turn  [6],  the  chief  orbit  is  most  often  not  circular.  A  dominant  downside  to  relaxing 
the  circular  chief  assumption  is  the  transition  from  a  time-invariant  system  to  a  time- 
varying  system,  in  that  a  spacecraft’s  dynamic  rates  are  explicit  functions  of  time. 
Realizing  the  future  impact,  in  the  decade  following  the  Clohessy- Wiltshire  publica- 


tion,  an  independent  effort  was  undertaken  by  Lawden  [14]  along  with  Tschauner  and 
Hempel  [15].  Tschauner  and  Hempel  were  pursuing  the  closed-form  solution  for  an 
arbitrary  elliptic  chief,  while  Lawden  was  attempting  to  describe  the  primer  vector 
with  respect  to  optimal  trajectories.  The  two  models  proved  to  be  nearly  identical, 
but  cumbersome  with  the  inclusion  of  the  below  integral  [12] 


I 


de 

sin^  9{1  +  e  cos  Oy 


(2.1) 


where  6  is  the  true  anomaly  of  the  chief.  The  integral  in  Eq.(2.1)  is  obviously  singular 
for  true  anomalies  of  integer  tt  multiples.  Carter  later  transformed  this  integral  to  the 
form  [16] 
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which  remains  nonsingular.  More  recently,  in  2002  Yamanaka  and  Ankersen  [6]  used 
basic  orbital  mechanics  relations  to  provide  a  solution  to  this  integral  as  a  function 
of  orbital  angular  momentum.  The  Yamanaka- Ankersen  model  will  be  discussed  in 
more  depth  further  in  Section  3.1.2. 


Alfriend  and  Yan  contend  [17]  that  for  the  nearly  thirty  years  following  the 
Tschauner- Hempel  and  Lawden  derivations,  little  work  was  done  in  the  held.  Some 
work  such  as  Carter’s  reformulation  of  the  singular  integral  [16]  was  performed,  but 
it  was  not  until  the  concept  of  satellite  formation  hight  was  popularized  in  the  late 
1990s  and  early  2000s  that  the  quest  for  model  hdelity  was  pursued.  The  quest  for 
high  hdelity  in  the  model  is  a  direct  consequence  of  desired  long  term  accuracy.  Most 
often,  the  unperturbed  Keplerian  orbit  assumption  is  made  to  soften  the  mathemat¬ 
ics.  However,  for  low  Earth  orbits  (LEO),  air  drag  and  oblateness  impact  the  orbit 
signihcantly  (See  Sections  2.5.1  and  2.5.2  for  specihcs  on  the  impact  on  relative  mo¬ 
tion). 


A  dominating  perturbation  to  consider  for  relative  motion  is  the  J2  ehect.  Mak¬ 
ing  the  assumption  of  a  circular  chief,  Schweighart  and  Sedwick  [7]  incorporated  an 
orbit  averaged  J2  modihcation  to  the  classic  HCW  model  that  captures  the  ehects  very 
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well.  The  HCW  equations  have  also  been  further  modihed  to  include  both  linear  [18] 
and  quadratic  [19]  drag  effects.  Chen  and  Jing  have  recently  developed  differential 
equations  including  both  air  drag  and  J2  on  the  arbitrary  elliptical  chief  [20]  using 
Lagrangian  dynamics,  but  has  no  closed  form  solution. 

Modeling  the  elliptical  chief  is  significantly  more  challenging  than  the  circular, 
reinforced  by  the  fact  that  the  Chen-Jing  model  does  not  have  a  closed  form  solution. 
In  recent  literature,  the  concept  of  orbital  element  differences  to  describe  the  relative 
motion  has  come  to  light.  Rather  than  expressing  the  trajectory  in  Cartesian  coordi¬ 
nates  {x,y,z),  the  more  intuitive  idea  of  using  the  difference  in  the  chief  and  deputy 
classical  or  non-singular  orbital  elements  can  be  used  to  propagate  to  a  future  state. 
Alfriend  has  shown  that  relative  motion  theories  using  these  states  tends  to  higher 
accuracy  than  using  a  Cartesian  or  curvilinear  frame  [21].  In  this  manner,  the  com¬ 
plexity  of  the  model  increases  significantly  if  the  instantaneous  orbital  elements  are 
used,  and  is  simplified  if  the  model  operates  in  mean  orbital  element  space.  Schaub 
makes  use  of  this  in  [22]  with  the  propagation  existing  in  mean  orbital  element  space. 
Hamel  and  Lafontaine  [23]  derive  a  time-varying  state  space  form  using  these  differ¬ 
ences,  but  rely  on  an  estimated  time  of  flight  solution.  Gim  and  Alfriend  [8]  derive 
a  state  transition  matrix  that  progates  either  the  initial  Cartesian  conditions  or  the 
initial  orbital  element  differences  to  a  desired  time  using  either  the  mean  or  osculating 
J2  effect.  The  GA  model  is  of  significant  importance  to  the  analytical  investigation 
in  this  study  and  is  described  in  depth  in  Section  3.1.1. 

There  do  exist  approaches  to  mapping  the  relative  motion  other  than  solving 
ordinary  differential  equations.  For  example,  Wiesel  applies  Floquet  theory  in  [5]  to 
the  relative  motion  problem.  In  this  work,  Wiesel  incorporates  all  zonal  harmonics  to 
the  HCW  model  and  produces  a  model  with  two  modes  that  are  linear  in  time,  and 
completley  describes  the  solution  by  a  periodic  vector  and  a  periodic  modal  matrix. 
As  another  example,  work  done  by  Kolemen  and  Kasdin  [24]  introduces  eccentricity  as 
a  perturbation  to  the  linearized  HCW  model  using  Hamiltonian  mechanics.  Another 
non-traditional  solution  is  that  of  Vadali  [25]  and  Sengupta,  Vadali,  and  Alfriend  [26] 
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in  which  the  motion  of  both  the  chief  and  deputy  are  normalized  by  their  respective 
radii  and  projected  onto  a  unit  sphere.  Following  the  projection,  spherical  trigonome¬ 
try  reveals  an  exact  kinematic  relation,  and  is  found  to  be  accurate  for  relative  ranges 
of  up  to  approximately  160  km  [17]. 

In  a  most  recent  effort  at  the  Air  Force  Institute  of  Technology  (AFIT),  Kirk 
Johnson  [27]  applied  a  ’’Virtual  Chief’  method  to  the  relative  motion  problem.  A 
circular  virtual  chief  on  an  orbit  known  a  priori  whose  orbital  elements  are  identical 
to  those  of  the  chief  with  the  exception  of  the  zero  eccentricity  is  propagated.  The 
motion  of  the  physical  chief  and  deputy  are  propagated  in  accordance  with  HCW 
dynamics  with  respect  to  the  virtual  reference.  The  coupling  of  linearization  errors  of 
both  the  chief  and  deputy  with  respect  to  the  HCW  propagation  introduced  significant 
error  when  compared  to  two  body  dynamics.  However,  the  possibility  of  retaining 
higher-order  terms  in  the  set-up  of  the  model  may  increase  the  accuracy  of  the  model. 
Johnson  also  proposed  in  his  thesis  a  set  of  six  parameters  that  defined  the  motion  of 
the  VC  trajectory  with  respect  to  geometry,  drift,  phasing,  periodicity,  and  skewness. 
As  a  side  study  to  this  thesis,  the  conversion  of  these  parameters  to  ROEs  under  a 
circular  assumption  is  presented  in  Appendix  E. 

Linearizations  made  also  impact  the  physical  space  of  the  model;  for  example, 
an  initial  condition  on  relative  displacement  for  one  model  does  not  necessarily  cor¬ 
respond  to  an  initial  condition  for  numerical  integration.  Therefore,  initial  condition 
matching  has  also  populated  the  literature  surrounding  this  topic.  Of  particular  in¬ 
terest  is  the  concept  of  creating  a  bounded,  periodic  relative  orbit  by  specifying  initial 
conditions.  The  problem  resulting  from  this  is  an  initial  condition  in  one  linearized 
model  does  not  necessarily  correspond  to  the  same  condition  in  another  linearized 
model.  Inalhan,  Tillerson,  and  How  present  an  algorithm  to  initialize  the  velocity 
space  as  a  function  of  position  at  any  point  in  the  orbit  [28].  Gurfil  has  recently 
adapted  the  periodicity  problem  in  terms  of  energy  matching  [29],  which  allows  ve¬ 
locity  to  be  extracted  from  a  quadratic  expression.  In  Schaub’s  text  [30],  an  in-depth 
discussion  is  made  concerning  period  matching  and  J2  invariant  relative  orbits.  From 
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this  brief  survey,  the  overarching  method  of  initial  condition  matching  is  in  the  form 
of  maintaining  equal  periods  between  the  two  orbits.  In  some  cases  this  may  relax  to 
equal  semi-major  axes;  however,  in  the  presence  of  the  J2  disturbance,  short-period 
oscillatons  in  the  chief  and  deputy  semi-major  axes  will  modify  the  period.  Moreover, 
in  the  presence  of  atmospheric  drag,  the  decay  of  the  orbit’s  energy  will  also  challenge 
energy  matching  between  the  chief  and  the  deputy. 

2. 2  Coordinate  Frames 

Now  being  familiar  with  the  importance  and  applicability  of  relative  satellite 
motion,  it  is  now  necessary  to  introduce  the  environment  which  the  majority  of  models 
describe. 

2.2. 1  Inertial  Reference  Frame.  A  geocentric  reference  frame  will  be  used  for 
the  purpose  of  the  current  study.  The  Earth-centered-inertial  frame  originating  from 
the  center  of  the  Earth  is  the  frame  of  choice.  Seen  in  Fig.  2.1,  the  I  vector  is  in  the 
direction  of  the  vernal  equinox,  the  K  vector  is  vertical  through  the  north  pole,  and 
the  J  vector  is  orthogonal  eastward  to  I  and  lies  in  the  Equatorial  plane.  Serving  as 
one  of  the  most  common  frames  in  orbital  mechanics,  it  is  deemed  “inertial  enough.” 
Slight  changes  occur  in  the  location  of  the  vernal  equniox  and  in  the  equitorial  plane; 
however,  speficiation  of  a  certain  epoch  often  implies  an  inertial  system. 

2.2.2  Local  Vertical- Local  Horizontal.  The  majority  of  the  relative  motion 
models  studied  express  their  solution  in  the  local-vertical,  local-horizontal  (LVLH) 
frame.  This  is  a  frame  that  assumes  the  satellite  is  a  point  mass  and  rotates  at  the 
same  angular  rate  as  the  satellite  in  its  orbit.  For  example,  a  LVLH  frame  attached 
to  an  unperturbed  circular  chief  will  rotate  exactly  at  the  mean  motion,  while  a 
LVLH  frame  attached  to  an  unperturbed  non-circular  chief  will  rotate  at  the  time 
derivative  of  the  true  anomaly.  Another  term  for  this  frame,  and  used  interchangably 
throughout  this  study,  is  the  RIC  (radial,  in-track,  cross-track)  frame.  The  radial 
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Figure  2.1:  Earth  Centered  Inertial  Reference  Frame 

(6^)  axis  is  in  the  direction  of  the  position  vector  from  the  center  of  the  Earth  to  the 
satellite,  positive  outward.  The  cross-track  (oh)  axis  is  parallel  to  the  orbital  angular 
momentum  vector  in  the  orbit  normal  direction.  The  in-track  (og)  axis  completes  the 
orthogonal  set.  Equation  2.3  expresses  the  LVLH  vectors  mathematically,  following 
Schaub  [22]. 

dr  =  r/r 

dh  =  h/h  (2.3) 

Off  =  Oh  X  dr 

where  r  is  the  position  vector  from  the  center  of  the  central  body  to  the  chief,  h  is 
the  specihc  angular  momentum  vector  (h  =  r  x  r),  and  the  unbolded  scalars  are  the 
2-norm  of  the  respective  vectors.  Figure  2.2  provides  a  visualization  of  the  reference 
frame 

The  location  of  the  deputy  in  the  relative  LVLH  frame  must  now  be  discussed. 
The  chief,  although  constantly  moving  and  rotating,  is  kept  at  the  origin  of  the 
frame.  The  position  vector  from  the  chief  to  the  deputy  is  denoted  as  p,  and  its 
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of  inertial  position  vector 

of  angular  momentum 

Og  completes  orthogonal  triad;  spans  in¬ 
plane  motion  with 

Inertial  Oihit 

Figure  2.2:  Local  Vertical  Local  Horizontal  Reference  Frame 

time  derivative  as  p.  The  Cartesian  components  of  the  relative  position  and  velocity 
vectors  are  expressed  in  the  radial  (x),  in-track  [y),  and  cross-track  [z)  directions.  For 
this  study,  the  following  notation  in  Eq.  2.4  is  used  for  the  relative  position  vector 

^  X 

P=  y 

V 

The  relative  velocity  vector  is  presented  similarly,  only  using  scalar  relative  velocity 
components  where  the  relative  position  states  are.  Figure  2.3  provides  a  visualization 
of  the  deputy’s  location  relative  to  the  chief.  It  is  important  now  to  note  that  even 
though  a  time  history  of  a  relative  orbit  will  show  movement  around  the  origin,  the 
deputy  is  in  reality  not  orbiting  around  the  chief.  The  central  body  remains  the  source 
of  the  primary  gravitational  force  held. 

2.2.3  Curvilinear  Reference  Frame.  Similar  to  the  LVLH  frame  described  in 
2.2.2,  the  curvilinear  frame  is  nearly  identical  in  form.  The  Gim-Alfriend  model  uses 
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Figure  2.3:  Relative  Position  Vector  in  the  LVLH  Frame 

curvilinear  coordinates  to  obtain  more  accurate  results  [8] .  The  calculated  states  still 
represent  the  same  physical  phenomenon,  and  may  still  be  used  for  ROE  substitution. 
The  X  coordinate  is  taken  as  the  difference  in  radii  between  the  chief  and  the  deputy, 
the  y  and  2:  coordinates  are  taken  as  curvilinear  distances  along  and  perpendicular 
to  an  instantaneous  imaginary  circle  on  the  reference  orbital  plane.  In  this  sense,  the 
relative  position  vector  can  be  thonght  to  “bend”  along  the  trajectory. 

2.3  The  Clohessy- Wiltshire  Model 

The  most  famous  relative  motion  model  is  that  of  Clohessy  and  Wiltshire  [4].  As 
the  HCW  model  is  the  basis  for  the  relative  orbit  element  realization,  the  following 
section  presents  a  beginning  to  end  derivation  and  trend  analysis  as  discnssed  in 
Vallado  [31]. 

2.3.1  Derivation.  To  begin  the  derivation  of  the  HCW  model,  a  problem 
statement  is  hrst  defined.  It  is  desired  to  determine  the  analytical  solution  for  the 
relative  motion  of  the  deputy  with  respect  to  a  circular  chief.  The  inertial  chief  and 
deputy  position  vectors  are  denoted  as  fy  and  fd,  respectively.  The  inertial  relative 
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range  vector  (*p)  is  now 


(2.5) 


■Tc 

If  one  desires  to  obtain  the  eqnations  of  motion,  the  inertial  relative  acceleration 
vector  is  needed  (Eq.  2.6) 

"p  =  rd-  fc  (2.6) 

Assnming  unperturbed  Keplerian  two-body  motion,  the  individual  accelerations  are 
well  known  and  can  be  represented  as  Eq.  2.7 


r  = 


(2.7) 


Now,  substituting  Eq.  2.7  in  Eq.  2.6,  the  differential  equation  becomes 


P  = 


pr-d  pvc 
'  ^3  ^  ^3 


r 


d 


(2.8) 


Having  the  basic  set-up,  there  is  now  a  need  to  simplify  the  system.  Using  the  law  of 
cosines,  the  magnitude  of  the  deputy’s  inertial  position  vector  is  then  found  as 


r 


2 

d 


J’c  +  y  +  2r). p 


(2.9) 


Substituting  Eq.  2.9  into  a  form  similar  to  the  two-body  acceleration 


U  _  U  +'  p 

'f'd  {rl  +  p‘^  +  2fc  •*  p)  t 


(2.10) 


Now  is  where  the  second  assumption  (the  hrst  being  unperturbed  Keplerian  motion) 
is  made.  It  is  assumed  that  the  relative  position  vector  is  much  smaller  than  the 
position  vector  of  the  chief.  Therefore,  p^  «  r^.  Applying  this  assumption  to  Eq. 
2.10  and  simplifying,  Eq.  2.11  falls  out 
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(2.11) 
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A  binomial  expansion  of  the  form 


{l  +  xY  =  ^{  =  l  +  nx+  +  ... 

k=0  ^  ' 


(2.12) 


is  then  applied  to  the  denominator  of  Eq.  2.11,  and  truncated  at  the  first  term  to 
arrive  at 


2r-;-vy^/^  3(2r-;--p) 

rl  )  2  r2 

Substituting  Eq.  2.13  and  substituting  in  Eq.  2.11, 


(2.13) 


fd  _  n,  +V  q  _  5  (2'!!  P) 


2  r'i 


(2.14) 


Taking  Eq.  2.14  and  substituting  in  Eq.  2.8 


P  = 


prd  prc 
^  d 


=  -p 


Tc  +*  P 


d  3(2r-;.V) 
2  r? 


+ 


prc 


(2.15) 


Noting  opposite  signs  and  simplifying,  we  arrive  at 


P  j.?, 


H  [  3r;2r;.V 


2rc  Tc 


+  V 


(2.16) 


It  is  now  noted  that  the  inertial  position  vector  can  be  expressed  in  the  LVLH  frame 
from  the  dot  products  in  Eq.  2.16.  This  follows  mathematically  as 


3  fc  2rc  p 

2  Tr  Tr 


=  —3dr  {dr  p) 


(2.17) 


=  —3xdr 


After  substitution,  we  arrive  at 


P 


P=  — ^(-3a;6,.  +*  p) 


(2.18) 
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As  this  is  the  inertial  acceleration,  and  the  desired  solution  is  in  the  LVLH  frame,  the 
transport  theorem  is  now  applied.  For  completion,  the  transport  theorem  is 

•|(,)=‘|(.)+-“x(.)  (2.19) 

where  the  superscript  b  denotes  a  non-inertial  reference  frame,  and  the  superscript  i 
denoted  an  inertial  frame.  Because  this  is  an  acceleration,  the  transport  theorem  is 
applied  twice  to  yield,  among  others,  the  well  known  Coriolis  and  centripetal  terms. 
The  full  transport  theorem  is  expressed  as 

r  +  ci}  X  r  +  2(i;  X  r  +  cJ  X  (cJ  X  r)  (2.20) 

Assuming  a  circular  chief,  the  angular  velocity  of  the  LVLH  frame  with  respect  to  the 
inertial  frame  is  the  mean  motion  of  the  chief.  That  is,  u  =  ^  r^Oh  =  noh  and  cJ  =  0. 

y 

Applying  Eq.  2.20  and  the  angular  velocity  to  Eq.  2.18,  the  hnal  vector  equation  of 
motion  is  found  as 

p  =  —n^{p  —  3xdr)  +  2  {nydr  —  nxdg)  xdr  yde  (2.21) 

Expressing  the  p  vector  in  component  form 

xdr  +  yde  +  zdh  =  —n^{{x  —  3x)dr  +  yd0  +  zdh)  +  ‘2.{nydr  —  nxdo)+n‘^xdr  +  n‘^ydg  (2.22) 

And  finally  collecting  the  vector  components,  the  set  of  three  linear,  coupled,  ordinary 
differential  equations  is  found  as 


0  =  X  —  2ny  —  3n^x 

0  =  y  +  2nx  (2.23) 

0  =  ^  +  71^  Z 
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It  it  worth  noting  that  if  an  external  force  is  present,  the  left-hand  size  of  zeros  in 
Eq.  2.23  are  replaced  by  more  specihc  expressions.  Finally,  to  place  the  differential 
equations  in  a  more  comfortable  matrix  form,  the  system  above  can  be  expressed  in 
state-space  form  as 


d) 

0 

0 

0 

1 

0 

0 

(x\ 

y 

0 

0 

0 

0 

1 

0 

y 

z 

0 

0 

0 

0 

0 

1 

z 

X 

3n^ 

0 

0 

0 

2n 

0 

X 

y 

0 

0 

0 

—2n 

0 

0 

y 

0 

0 

— 

0 

0 

0 

vJ 

(2.24) 


Having  the  familiar  form  of  X  =  AX,  basic  concepts  of  linear  systems  can  be  used  to 
hnd  the  closed-form  solutions  in  Section  2.3.2. 


Without  replicating  the  derivation,  Schaub  [30]  presents  the  coupled,  non-linear 
system  of  differential  equations  as 


x-2f  (ij-  y—^  -xp  -  ^  =  -4(^c  +  x) 

\  rj  ri  rl 

y +  2f  (x  -  -yp  = -^y  (2.25) 

V  ^cj  n 

y 

Z  = - wZ 

where  /  is  the  true  anomaly  of  the  chief.  This  system  remains  valid  for  arbitrarily 
large  orbits  with  a  non-circular  chief.  However,  making  the  same  assumptions  as 
earlier  (circular  chief,  close  proximity,  first  order  truncation)  this  model  does  in  fact 
yield  the  same  linear  time-invariant  set  given  in  Eq.  2.24  (the  HCW  model). 

Finally,  we  can  change  the  variables  to  a  curvilinear  frame  via  the  following 


X  =  6r 


y  =  rc66 


(2.26) 
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where  6r  and  66  are  differential  displacements  in  the  radial  and  angular  directions. 
The  HCW  equations  can  be  expressed  as 

0  =  (5r  —  2nrc69  —  3n^6r 

0  =  TcSd  +  2n5r  (2.27) 

0  =  ^  +  71^  Z 

The  use  of  curvilinear  coordinates  is  often  found  to  be  more  accurate  as  the  system 
will  naturally  bend  the  in-track  vector  along  the  orbital  path  [30]. 

2.3.2  Solving  the  HCW  Differential  Equations.  Now  having  a  system  of 
linear,  time-invariant  ordinary  differential  equations  (LTI  ODEs),  fundamental  prop¬ 
erties  of  linear  systems  can  be  applied  to  hnd  the  solution  to  Eq.  2.24.  The  hrst 
fundamental  principle  (applied  without  proof)  is  that  the  solution  to  the  unforced 
system  of  LTI  ODEs 

i  =  Ax  (2.28) 

where  A  is  the  plant  matrix,  and  x  is  the  vector  of  states,  is  of  the  form 


x(t)  =  <I>(t  -  to)xo 


(2.29) 


where  <h(t  —  to)  is  termed  the  state  transition  matrix  (STM).  As  STMs  form  a  sig- 
nihcant  foundation  for  this  study,  a  brief  survey  is  provided  in  Appendix  A.  To  hnd 
<h(t  —  to),  the  matrix  exponential  is  used 


$(t  -  to)  = 

*h(t  -  to)  =  (si  -  A)~^ 


(2.30) 
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with  s  being  the  Laplace  variable.  Substituting  in  the  plant  matrix  from  Eq.  2.24, 
the  STM  is  expressed  as 


s  0  0 

0  s  0 

0  0s 

-3n2  0  0 

0  0  0 

0  0 


-10  0 
0-10 
0  0-1 
s  —2n  0 

—2n  s  0 
0  0s 


(2.31) 


where  £  ^  is  the  inverse  Laplace  transform.  Without  replicating  the  matrix  algebra, 
Eq.  2.31  reduces  to  a  very  elegant  solution  of  the  form 


-  to)  = 


4  —  3  cos  6 

0 

0 

sin  0 

n 

f(l  -COS0) 

0 

6(sin6'  —  9) 

1 

0 

f(cos0-  1) 

-  sin  6^  -  — 

n  n 

0 

0 

0 

cos  9 

0 

0 

sin  6 

n 

Snsind 

0 

0 

cos  9 

2sin0 

0 

6n(cos9  —  1) 

0 

0 

-2  sin  9 

4  cos  9  —  3 

0 

0 

0  - 

-nsin9 

0 

0 

cos  6* 

:  —  to)-  The  state 

at  any 

time  is  now 

a  linear  combination  < 

(2.32) 


initial  conditions.  Often  in  the  literature  (for  example,  [32]),  the  HCW  matrix  is 
partitioned  into  four  3x3  matrices,  as  the  following 


<h(t  -  to) 


^rr  ^rv 


(2.33) 


where  the  subscripts  r  and  v  correspond  to  position  and  velocity,  respectively.  Having 
a  complete  closed  form  solution  to  the  LTI  ODEs,  the  system  behavior  can  now  be 
investigated  in  Section  2.3.3. 
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2.3.3  System  Behavior.  In  this  section,  only  a  brief  survey  is  given  on  the 
system  behavior  of  the  HCW  equations.  This  is  a  direct  consequence  of  the  ROEs 
completely  describing  the  trajectory.  As  such,  to  avoid  duplicate  derivatons,  only 
top-level  system  qualities  are  detailed,  and  the  more  mathematical  derivations  are 
detailed  in  Section  2.4. 

The  most  elegant  fallout  is  the  geometry  of  the  in-plane  motion  of  the  deputy 
relative  to  the  chief.  It  will  later  be  shown  that  the  radial  and  in-track  trajectory  is  a 
2:1  ellipse.  That  is,  the  in-track  motion  oscillates  at  twice  the  magnitude  of  the  radial 
with  an  orthogonal  phasing.  The  ellipse  is  centered  at  a  constant  radial  displacement, 
but  an  in-track  displacement  that  may  drift  if  certain  initial  conditions  are  not  met. 
Examining  Eq.  2.32,  the  only  secular  terms  are  in  the  y(t)  equation  is 

Secular  Term  =  {—6nxo  —  3yo)  t  (2.34) 

and  to  eliminate  the  drift  in  the  in-track  direction,  Eq.  2.34  must  be  null,  implying 

yo  =  -2nxo  (2.35) 

which  will  produce  a  2:1  ellipse  centered  at  a  constant  radial  and  in-track  displace¬ 
ment. 

The  cross-track  [z)  motion  is  completely  uncoupled  from  the  radial  and  in-track 
motion  in  the  linearized  case.  The  result  of  this  uncoupling  is  that  any  cross-track 
displacement  will  simply  superimpose  an  oscillatory  motion  on  the  in-plane  trajectory 
in  the  cross-track  direction.  Vallado  [31]  summarizes  of  the  effect  of  initial  conditions 
on  the  motion,  detailed  below 

1.  The  deputy  will  progress  further  in-track  if  beginning  in  a  lower  orbit  (xq  <  0). 

For  unperturbed  motion,  a  lower  orbit  implies  a  higher-energy  orbit.  The  deputy 

will  move  at  a  greater  relative  velocity  than  the  chief  and  move  ahead  in-track 
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2.  The  deputy  will  regress  further  in-track  if  beginning  in  a  higher  orbit  (xq  >  0). 
For  unperturbed  motion,  a  higher  orbit  implies  a  lower-energy  orbit.  The  deputy 
will  move  at  a  lower  relative  velocity  than  the  chief  and  move  behind  in-track 

3.  As  the  initial  radial  displacement  is  increased,  the  relative  motion  trajectory 
will  increase.  As  shown  later  in  the  ROE  development  the  size  of  the  relative 
orbit  is  determined  by  the  radial  displacement,  and  the  velocities  in  the  radial 
and  in-track  direction. 

4.  The  ellipse  can  be  centered  at  a  constant  displacement  and  the  relative  orbit 
made  periodic  if  the  initial  condition  yo  =  —‘2nxQ  is  enforced. 

5.  For  a  hxed  radial  displacement,  any  in-track  displacement  will  produce  the  same 
relative  motion,  as  the  deputy  will  still  be  in  the  same  orbit  at  the  same  energy. 

6.  A  cross-track  displacement  will  simply  superimpose  oscillations  over  the  in-plane 
motion. 

Now,  having  a  solid  foundation  for  the  HCW  motion,  the  relative  orbit  elements 
can  now  be  derived  and  shown  as  a  practical  parameterization  of  the  model. 

2.4  The  Relative  Orbit  Element  Realization 

Taking  the  results  from  Section  2.3,  the  relative  orbit  elements  are  now  derived. 
Being  that  the  ROEs  are  the  basis  of  this  study,  great  care  is  taken  in  their  derivation 
following  the  results  of  Lovell  [33] . 

2.4-1  Derivation.  The  ROEs  are  a  set  of  six  parameters  that  describe  the 
relative  motion  in  the  LVLH  frame  for  the  circular  chief.  The  six  ROEs  are  the 
semi-major  axis  (oe),  the  radial  displacement  (x^),  the  in-track  displacement  (y^), 
the  in-plane  phasing  (/5),  the  maximum  cross-track  amplihcation  {zmax),  and  the  out- 
of-plane  phasing  {'ip).  The  derivation  will  be  divided  into  the  ROEs  describing  the 
in-plane  and  the  out-of-plane  motion. 
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2.4-1 -1  In-Plane  ROEs.  The  two  ROEs  most  easily  derived  are  the 
radial  and  in-track  displacements.  Examining  the  HCW  expressions  for  the  motion 
in  these  two  directions 


a;=  ^sin0-  faxoT— ^  cos0+  (4x0  +  —'] 

n  \  n  J  \  n  J 

y  =  ( 6x0  -\-  sin  9  -\-  cos  9  —  {6nxo  -\-  3yo)t  +  ( yo  —  (2.36) 

\  n  J  n  \  n  J 


it  remains  obvious  that  the  x  expression  is  oscillatory  with  the  exception  of  an  offset 
term,  and  that  the  y  expression  is  oscillatory  with  the  exception  of  an  offset  and  a 
secular  term.  These  offsets  are  dehned  as  the  radial  and  in-track  displacements  such 


that 


Xd=  (  4x0  +  — 
n 


yd  =  -(6nxo  +  3yo)t  +  [yo-  — 

n 


(2.37) 


Letting  the  constant  term  (yo  —  ^)  =  ydo-,  and  recognizing  that  the  expression  for  yd 
is  a  linear  combination  of  Xd,  it  is  then  found  that  the  displacements  can  be  written 
as 


/I  ^  ^^0 

Xd  =  4xo  H - 

n 

yd  =  ydo  -  ^nxd{t  -  to) 


(2.38) 


To  determine  the  magnitude  of  the  oscillatory  in-plane  motion,  one  can  apply  the 
Hamonic  Addition  Theorem  (or  the  linear  combination  of  trigonometric  functions), 
which  states  that  a  function  written  as 


f  =  Asm9  I-  B  cos 9 


(2.39) 
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can  be  written  in  the  following  two  manners 


/ 

/ 


sin 


1^0  +  tan  ^ 


—  VA^  +  B^  cos 


1^0  +  tan  ^ 


(2.40) 


Applying  Eq.  2.40  to  the  radial  motion  in  Eq.  2.36,  and  letting  A  =  ^  and  B  = 
—  (3x0  +  we  see  that 

+  (3:.„  cos  (^»  +  ton-‘  (2-«) 

For  simplicity,  the  constant  phasing  parameter  /3o  can  now  be  defined  as 

/3o  =  tan"^ - (2.42) 

(3X0  +  ^)  ^  ^ 

Rewriting  the  radial  motion  using  C  =  +  (3xo  +  such  that 

X  =  —C  cos{9  +  /do)  (2.43) 


The  argument  in  the  cosine  function  can  also  be  reduced  using  (3  =  6  +  (3o^  providing 
the  final  radial  expression 

x  =  — Ccos/d  (2.44) 

Having  a  simplihed  expression  for  the  radial  motion,  the  in-track  motion  can  be 
simplihed  in  the  same  manner.  Applying  Eq.  2.40  to  the  in-track  expression,  and 
letting  A  =  (6xo  +  ^)  and  B  =  (^),  it  follows  that 
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Substituting  the  phase  and  dummy  magnitudes, 


y  =  2Csm(3 


(2.46) 


Thus  it  follows  that  the  in-track  motion  oscillates  orthogonal  to  the  radial  motion 
at  twice  the  magnitude.  This  corresponds  to  the  dehnition  of  a  2:1  ellipse  with  the 
semi-major  axis  on  the  in-track  axis.  The  in-plane  ROEs  have  now  been  described 
and  are  summarized  below,  removing  the  subscripts  other  than  on  the  initial  ROEs: 


(Xq 


Xd  =  -\ - 

n 

3 

Vd  VdO 


P  =  6  +  tan 


-1 


X 


3nx  +  2y 


—  9  +  (3^ 


(2.47) 


where  9  =  n{t  —  to).  The  transformation  to  Cartesian  position  coordinates  was  also 
found,  while  the  velocity  coordinates  are  simply  the  time  derivatives 


X  =  cos  fJ  -F  Xd 

y  =  aeSin/3  -1-  yd 

Ctefl  . 

X  =  - sm  n 

2 

Stz 

y  =  tteU  cos  /3  -  —Xd 


(2.48) 


2.4-1 -2  Out-of-Plane  ROEs.  As  the  out-of-plane  motion  is  uncoupled 
from  the  in-plane,  only  two  parameters  (other  than  the  mean  motion  and  time)  are 
needed  to  fully  describe  the  motion.  Examining  the  cross-track  position 

z  =  zqcos9  ^ — ^sin6'  (2.49) 

n 
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it  is  fairly  obvious  the  motion  is  solely  sinusoidal  with  no  offset.  Directly  applying 
the  Harmonic  Addition  Theorem  to  Eq.  2.49,  and  letting  A  =  ^  and  B  =  zq,  the 
cross-track  motion  is  written  as 


z  =  \l  {  —  ^  -I-  2:n  sin  f  6*  -I-  tan  ^  f  ^ 
n 


£0 

n 


(2.50) 


Similar  to  the  in-plane  motion,  initial  and  current  phase  angles  ('^o;  t/’)  can  be  dehned 
such  that 


^0  =  tan-i ^ 

^  =  e  +  iJo 

and  the  cross  track  motion  reduces  to 


(2.51) 


2:  = 


—  -1-2:0  sin  Ip 


(2.52) 


The  motion  has  now  been  reduced  to  a  simple  sinusoid,  with  an  amplitude  of 


^m.ax  — 


-I 

n  ' 


(2.53) 


The  cross-track  motion  can  now  be  specified  as 


^  =  Zmax  sin  tp 


(2.54) 


Consequently,  the  two  ROEs  describing  the  cross-track  motion  can  be  described  as 


-I  +4 

n 


(2.55) 


ip  =  6  +  ipQ  =  6  +  tan  ^ 

V  i 
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where  6  =  n{t  —  to)-  The  inverse  transformation  to  Cartesian  coordinates  is 


2:  =  Zjnax  sin  ijj 
Z  =  ZraaxU  COS  tjj 


(2.56) 


2.4-2  System  Behavior.  The  time  evolution  of  the  relative  orbit  elements  is 
of  signihcant  importance  to  this  study.  The  linear  time  invariant  system  effected  by 
the  HCW  assumptions  yields  an  elegant  time  history  for  the  ROEs.  However,  when 
this  parameterization  is  applied  to  a  perturbed  environment  or  an  elliptical  chief, 
there  is  no  reason  to  expect  the  response  to  remain  the  same. 

Immediately  from  the  expressions  for  the  radial  and  in-track  motion,  the  ex¬ 
pression  for  Oe  remains  constant  over  time.  Also  from  the  radial  motion,  the  radial 
displacement  remains  constant  over  time.  The  in-track  motion  will  drift  at  a  rate  of 


Vd  = 


3n 

3n 


Xd 


2 


(2.57) 


which  can  now  be  utilized  for  a  boundedness  condition.  If  Xdo  =  0,  Eq.  2.57  reduces 
to  zero,  and  the  relative  ellipse  will  remain  stationary.  The  phase  angle  (3  increases 
linearly  with  time  from  the  initial  state  j3o  by  an  amount  proportional  to  the  mean 
motion  of  the  chief.  The  cross-track  ROEs  follow  a  similar  analysis.  The  maximum 
cross-track  amplihcation  remains  at  a  hxed  value,  and  the  cross-track  phasing  behaves 
similarly  to  the  in-plane  phasing. 

The  signihcance  of  the  ROE  parameterization  is  that  the  relative  motion  of  the 
deputy  can  be  found  as  an  instantaneous  elliptical  path  that  is  centered  at  yd,  0) 
[33].  If  the  relative  orbit  is  drifting,  it  will  drift  linearly  as  specihed  by  the  expression 
for  yd-  The  phase  angle  f3  is  related  to  the  angular  position  to  the  deputy  in  the 
relative  orbit.  The  cross-track  motion  intersects  the  LVLH  frame  of  the  chief  when 


the  phasing  is  an  integer  multiple  of  tt,  and  is  takes  on  the  value  of  Zmax  and  —Zmax 
when  the  sine  function  is  maximized  {ip  =  7r/2,37r/2,  respectively). 

If  the  angle  'y  =  ip  —  j3  is  used  as  the  constant  phase  difference  between  the 
cross-track  and  in-plane  motion,  then  the  cross-track  expressions  can  be  written  as 


2:  =  Zmax  sin(7  +  P) 
z  =  Zmaxn  cosi'j  +  P) 


(2.58) 


and  the  only  time-varying  ROEs  are  P  and  yd-,  both  of  which  change  linearly  with 
time. 


2.4.3  Summary.  To  summarize  the  ROEs,  analytical  expressions  for  the 
time  history  is  provided  as  detailed  in  Lovell  [33],  as  well  as  the  Cartesian  conversion. 
The  ROEs  can  be  coverted  from  Cartesian  coordinates  as 


(Xg 


Xd 


Ud 


Zmax 


Ip 


tan-  (-) 


(2.59) 
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The  set  R  is  now  defined  for  future  use  as  R  =  [og,  /?,  t/’] 

evolution  of  the  ROEs  is 

dg  cigo 

3/?/ 

Dd  DdO  ~^^di 

p  =  Po  +  nt 

^max  —  ^maxO 

Ip  =  ipQ  +  nt 

The  Cartesian  conversion  is 

a;  =  —  —  cos  p  +  Xd 
y  =  aeSin/3  +  yd 
Z  =  Zmax  sin  p: 

.  n 

X  =  — nsinp 
2 

y  =  tteU  cos  p  —  -nxd 
z  =  Zmaxn  cos  pj 

A  visualization  of  the  the  in-plane  ROEs  is  provided  in  Fig.  2.4  with  (3*  serving 
as  the  physical  interpretation  of  the  phase  angle. 

A  visualization  of  the  out-of-plane  ROEs  is  provided  in  Fig.  2.5  with  7*  serving 
as  the  in-plane  physical  interpretation  of  the  relative  ascending  node. 

2.5  Theoretical  Impacts  of  Eccentricity  and  J2  on  the  Relative  Trajec¬ 
tory 

Before  delving  into  the  analytical  and  numeric  results,  the  following  sections 
present  the  actual  impacts  expected  in  the  analytics.  The  effect  of  chief  eccentricity  is 
explored.  The  J2  perturbation  is  initially  surveyed,  and  its  application  as  a  perturbing 


The  time 


(2.60) 


(2.61) 
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a 


In-track 


Figure  2.4:  Visualization  of  the  In-Plane  Relative  Orbit  Using  ROEs 


Cross-track 


Figure  2.5:  Visualization  of  the  Out-of-Plane  Relative  Orbit  Using  ROEs 
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force  on  relative  motion  detailed.  It  should  be  noted  that  the  combination  of  both 
eccentricity  and  J2  will  induce  non-linear  combinations  of  the  impacts  described  in 
this  section. 

2.5.1  The  Impact  of  a  Non-Circlar  Chief  Orbit  on  Relative  Motion.  Jiang, 
et.  ah  [34]  investigated  the  unperturbed  hrst-order  relative  motion  to  a  signihcant 
degree  using  algebraic  methods.  It  was  concluded  that  the  radial/in-track  projection 
of  the  orbit  may  have  no  more  than  one  self-intersection  per  orbit,  and  that  the 
in-track/cross-track  and  radial/cross-track  may  interesct  no  more  than  two  times 
per  orbit.  Another  signihcant  conclusion  made  is  that  the  general  three-dimensional 
motion  is  only  planar  for  degenerate,  unrealistic  cases.  The  relative  motion  was 
found  to  rest  on  three  quadric  surfaces:  a  hyperboloid  of  one-sheet,  an  elliptic  cone, 
or  an  elliptic  cylinder.  Sengupta  [35]  utilizes  the  solution  to  the  Tschauner  and 
Hempel  equations  and  expands  the  in-track  and  cross-track  components  as  a  Fourier 
series  to  examine  the  effects  of  eccentricity.  Using  both  true  anomaly  and  time  as 
independent  variables,  Sengupta  expands  the  motion  using  Cauchy  residue  theory 
to  be  functions  of  constants,  a  dominant  harmonic,  size  parameters,  and  high-order 
harmonics.  Sengupta  notes  hve  signihcant  eccentricity  ehects  due  to  the  expansion 
using  true  anomaly: 

1.  The  presence  of  high-order  harmonics  in  the  motion  is  the  primary  mode  of  the 
deviation  from  the  2x1  trajectory. 

2.  The  amplitude  of  the  primary  harmonic  is  scaled  by  the  eccentricity,  which  re¬ 
sults  in  a  cross-track  expansion  and  in-track  decrease  as  eccentricity  is  increased. 

3.  A  phase  shift  is  induced  to  the  sinusoids  of  the  HCW  including  terms  of  (e^,  e®, ...). 

Ignoring  the  higher  order  eccentricity  terms  yields  the  original  phase. 

4.  Constant  terms  displace  the  center  of  the  relative  orbit,  and  is  dependent  on 
the  initial  in-plane  and  out-of-plane  phasing. 
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5.  The  relative  orbit  plane  is  skewed  by  the  higher  order  harmonics.  Most  notably, 
the  plane  of  symmetry  of  the  cross-track  motion  that  is  normal  to  the  in-track 
trajectory  is  skewed  as  a  function  of  the  eccentricity. 

Sengupta’s  expansion  in  time  exhibits  the  same  quantative  effects  as  the  true  anomaly 
expansion. 

2.5.2  The  Impact  of  the  J2  Perturbation. 


2.5.2. 1  Background  on  the  J2  Perturbation.  Motion  about  a  spherical 
central  body  dictates  a  constant  attraction  at  a  given  range  with  no  regard  to  angular 
location.  This  spherical,  point-mass  assumption  leads  to  an  inverse  square  gravity 
held.  Oblateness  effects  such  as  the  equatorial  will  dramatically  modify  the  gravi¬ 
tational  attraction,  and  for  an  arbitrary  body  the  differential  specihc  gravitational 
potential  acting  on  a  point  A  can  be  expressed  as  [30] 


dV 


Gdm 

s 


(2.62) 


with  G  as  the  universal  gravitational  constant,  dm  as  a  differential  mass  on  the  body, 
and  s  is  the  magnitude  of  the  position  vector  between  the  differential  mass  and  A. 
Following  the  introduction  of  Legendre  polynomials,  it  is  well  documented  that  the 
gravitational  potential  held  on  a  body  can  be  expressed  as  (  [30],  [31]) 


k=l 


Pk{cos'y)dm 


(2.63) 


where  r  is  the  norm  of  the  position  vector  from  the  body  centroid  to  the  point  A,  p  is 
the  norm  of  the  position  vector  from  the  body  centroid  to  the  dm,  7  is  angle  between 
p  and  r,  and  P^  is  the  Legendre  polynomial  such  that 


Pn+l{,J^) 


2n  +  l 
n  +  l 


uPn{ir) 


n 

n  -I-  1 


T’n-i(z^) 


(2.64) 


33 


(2.65) 


The  Legendre  polynomials  have  zero  mean  and  are  orthogonal  to  one  another. 

A  change  of  coordinates  from  Cartesian  to  spherical  will  result  in  an  inhnite  se¬ 
ries  that  expresses  the  gravitational  potential  in  terms  of  spherical  harmonics.  Rather 
than  expressing  it  as  Eq.  2.63,  the  gravitational  potential  can  eventually  be  written 
as 

h 

\  k=2 

where  Jk  is  the  zonal  gravitational  harmonic,  r^g  is  the  equatorial  radius,  and  0  locates 
the  position  vector  vertically  from  the  equatorial  plane.  Determination  of  the  zonal 
harmonics  is  typically  done  via  an  estimation  algorithm  extracted  from  observations. 
The  value  for  the  J2  harmonic  has  a  value  of  1082.63x10“®  [30]. 

Now,  using  the  classical  orbital  elements  a,  e,  i,  D,  u,  Mq,  numerous  partial  deriva¬ 
tives  and  extensive  use  of  Lagrangian  mechanics  will  yield  Lagrange’s  planetary  equa¬ 
tions,  presented  in  Eq.  2.66  using  Schaub’s  [30]  notation  but  placed  in  matrix  form. 
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where  a  and  b  are  the  semi-major  and  semi-minor  axes,  and  R  is  the  the  disturbance 
potential,  which  is  now  known  from  the  spherical  harmonic  gravitational  potential  in 
Eq.  2.65.  Truncating  the  disturbance  potential  at  the  dominant  second  harmonic,  R 
becomes 


J2fi  /re  ^2 


=  (3sinV-l) 

2  r  \  r 


(2.67) 


A  change  of  variables  (sin  0  =  sin(a;  -t-  z/)  sin  i)  then  allows  to  disturbance  potential  to 
be  written  as 

Ja/i  /re  ^  2 


R{r)  = - -  (—')  (3(sin(c<;  -|-  z/)  sini)^  —  1) 

2  r  V  r  / 
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Taking  the  gradient  of  Eq.  2.68,  the  osculating  time  rates  of  the  orbital  elements  can 
be  found  from  the  matrix  multiplication  in  Eq.  2.66.  The  J2  perturbation  induces 
both  short  and  long  term  periodic  effects,  in  addition  to  secular  drift  rates.  For 
relative  motion,  it  is  typically  assumed  that  the  differences  in  the  short  and  long  term 
variations  in  the  chief  and  deputy  are  negligible.  Thus,  it  is  the  secular  time  rate 
that  is  the  primary  perturbation  in  relative  motion.  An  orbit-average  of  Lagrange’s 
planetary  equations  using  the  J2  disturbance  yields  the  following  mean  orbital  element 
drift  rates  in  Eq.  (2.69  [31]. 


h  = 
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where  the  barred  quantities  indicate  mean  elements. 

The  resulting  orbit  averaged  COEs  make  intuitive  sense.  Using  conservation  of 
energy  arguments,  the  semi-major  axis  can  be  shown  to  vary  only  periodically  [31]. 
Inherent  in  the  construction  of  spherical  harmonics  is  the  assumption  of  an  axially 
symmetric  central  body;  this  implies  that  the  angular  momentum  about  the  polar 
axis  is  a  constant.  If  the  inclination  varied  secularly  from  J2,  at  one  point  this  term 
would  equate  to  zero,  which  also  directly  implies  that  the  eccentricity  can  only  vary 
periodically. 


2. 5. 2. 2  Effect  on  Relative  Motion.  Short  period  effects  of  J2  induce 
small  oscillations  in  the  orbit  elements,  but  do  not  contribute  to  the  orbital  drift  [30], 
as  errors  from  short  period  can  be  expected  on  the  order  of  50-100  m  [31]  and  nearly 
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identical  between  the  chief  and  deputy.  The  long  period  effect  is  the  rotation  of 
apsides  and  will  have  a  near  negligible  difference  between  the  chief  and  the  deputy 
for  close  proximity  operations. 

To  better  understand  what  effects  will  occur  from  the  perturbations,  the  linear 
transformation  from  orbital  elements  to  the  Hill  frame  (presented  in  [30])  is  examined. 
The  radial  displacement  is  merely  the  difference  in  the  two  orbital  radii.  From  the 
two-body  trajectory  relation,  the  orbital  radii  are  functions  of  a,  e,  and  true  anomaly. 
Using  mean  orbital  elements,  the  mean  change  in  the  radial  direction  will  be  unaffected 
by  a  and  e  oscillations,  but  will  be  impacted  by  the  mean  anomaly  difference  (which 
obviously  corresponds  to  a  true  anomaly  difference)  [30].  The  in-track  and  cross¬ 
track  motions  are  functions  of  the  chief’s  orbital  radius,  inclination,  and  true  latitude, 
and  also  are  functions  of  the  differences  in  the  true  latitude,  inclination,  and  right 
ascenscion  of  the  ascending  node.  Greater  variations  are  to  be  expected  in  these 
components  as  the  mean  rates  impact  this  motion  signihcantly. 

Schweighart  [36]  provides  a  well-rounded  discussion  concerning  where  the  J2 
effect  would  manifest  itself  in  the  relative  trajectory,  but  it  is  primarily  for  circular 
reference  orbits.  The  variance  in  the  ascending  node  modihes  the  amplitude  of  the 
cross-track  motion.  Schweighart  compares  the  cross-track  relation  to  a  ’’scissoring 
effect”  [36]  in  that  as  the  scissors  are  opened,  the  intersection  of  the  blade  tips  quickly 
moves  towards  the  handle,  but  opening  further  decreases  the  rate  of  the  intersection. 
The  blades  correspond  to  the  orbital  planes  and  their  movement  due  to  secular  drift  in 
the  ascending  node.  The  rotation  of  the  apsides  moves  the  location  of  closest  approach, 
modifying  the  period  and  minima  locations  of  the  radial  and  in-track  components. 
This  change  impacts  the  periodicity  of  the  radial  and  in-track  motion  in  addition  to 
modifying  the  period  of  the  reference  orbit.  A  similar  effect  is  found  from  the  secular 
change  of  the  mean  anomaly  which  also  modihes  the  period  of  the  reference  orbit. 

The  orbital  element  rate  of  the  mean  anomaly  and  argument  of  perigee  can  be 
collapsed  to  a  drift  in  the  mean  argument  of  latitude  {6  =  M  +  u).  The  difference 
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between  the  two  6  values  will  result  in  in-plane  angular  deviation.  The  differences 
in  will  lead  to  cross-track  amplihcation.  To  maintain  a  static  relative  orbit,  initial 
conditions  can  be  made  to  construct  invariant  orbits  based  on  equating  the  secular 
drift  of  the  latitude  and  ascending  node. 
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III.  Description  of  Applied  Models  and  Methodology 

As  a  means  to  provide  clarity  on  the  state  transition  matrices  and  models  applied 
throughout  the  remainder  of  this  study,  the  following  sections  will  present  in  detail 
the  methodology  and  results  of  the  applied  models.  The  Gim-Alfriend,  Yamanaka- 
Ankersen,  Schweighart-Sedwick,  and  numerically  integrated  truth  model  will  be  dis¬ 
cussed  in  detail.  The  application  of  each  model  to  this  study  is  then  discussed.  The 
processes  used  to  perform  the  operational  mapping,  stationary  orbit  initialization, 
and  the  guidance  algorithm  are  also  discussed. 

3. 1  Applied  Models 

The  following  section  details  the  propagation  in  the  Gim-Alfriend,  Yamanaka- 
Ankersen,  and  Schweighart-Sedwick  relative  motion  models.  The  Gim-Afriend  and 
Yamanaka-Ankersen  systems  are  presented  in  closed  form,  while  the  Schweighart- 
Sedwick  system  is  given  initially  as  a  set  of  differential  equations  that  are  later  solved 
by  the  author. 

3.1.1  Gim-Alfriend.  The  Gim-Alfriend  model  [8]  is  a  geometric  method 
to  determine  the  relative  motion  dynamics  to  include  both  chief  eccentricity  and 
the  J2  effect.  The  model  uses  the  non-singular  orbital  elements  defined  as  the  set 
e  =  (a,  0,  i,  gi,  g2;  where  a  is  the  semi-major  axis,  i  is  the  inclination,  9  is  the 
argument  of  latitude,  and  G  is  the  right  ascenscion  of  the  ascending  node.  The 
parameters  qi  and  q2  are  dehned  as 


gi  =  e  cos  00 
q2  =  e  sin  u 


(3.1) 


where  e  is  the  eccentricity,  and  00  is  the  argument  of  periapsis.  Using  the  argument 
of  true  latitude  and  the  g(.)  parameters  avoids  singularities  as  true  anomaly  and  the 
argument  of  perigee  are  undefined  for  circular  orbits.  For  higher  accuracy,  the  use  of 
curvilinear  coordinates  is  used  as  discussed  in  Section  2.2.3. 
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The  GA  model  provides  a  STM  to  propagate  the  Cartesian  initial  conditions 
forward  using  either  osculating  or  mean  orbital  elements.  Using  the  GA  verbiage,  the 
propagation  goes  as 

X(t)  =  [A{t)  +  aB{t)]6e{t)  (3.2) 

where  A  is  the  unperturbed  transformation  from  orbital  elements  to  Cartesian  coordi¬ 
nates,  B  is  transformation  from  orbital  elements  to  Cartesian  coordinates  containing 
J2  terms,  a  is  equal  to  3J2RI,  Re  is  the  radius  of  the  Earth,  and  Se(t)  is  the  instan¬ 
taneous  orbital  element  differences  ((5e  =  —  Oc).  The  vector  Se(t)  is  propagated  as 

(5e(t)  =  0e(t,to)5e(to)  (3.3) 

where  (pe  is  the  state  transition  matrix  for  orbital  element  differences.  The  propagation 
from  Cartesian  initial  conditions  to  the  current  Cartesian  state  is  thus 

X(t)  =  $j,(t,to)X(to)  ^  ^ 

(3.4) 

<hj2  =  [A{t)  +  aB{t)]4>e{t,to)[A{to)  +  aB{to)]~^ 

This  is  the  point  at  which  the  model  splits  between  the  use  of  osculating  orbital 
elements  and  mean  orbital  elements.  If  mean  orbital  elements  are  used,  the  <hj2  state 
transition  matrix  is  expressed  as  <hj2  such  that 

X(t)  =  l>j2(^Co)X(to)  ^  ^ 

(3-5) 

<hj2  =  [A{t)  +  aB{t)](j)e{t,to)[A{to)  +  aB{to)]  ^ 

In  Eq.  3.5,  A  is  equivalent  to  A;  however,  osculating  elements  invoke  a  different 
angular  velocity  and  B  ^  B. 

If  osculating  elements  are  to  be  used,  the  propagation  goes  as 
X(t)  =  4>j2(t,to)X(to) 

(3.6) 

=  [A{t)  +  aB{t)]D{t)(j)e{t,to)D  ^(to)[A(to)  +  a5(to)]  ^ 
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where  the  D  matrix  is  the  Jacobian  of  the  osculating  orbital  elements  with  respect  to 
the  mean  orbital  elements.  That  is 

m  =  (3.7) 

^^mean 

which  includes  the  short  and  long  period  oscillations  in  the  osculating  classical  orbital 
elements.  The  members  of  these  matrices  are  entirely  too  cumbersome  to  list  here 
and  are  detailed  in  the  appendices  of  [8]. 

Assumptions  inherent  in  the  GA  model  are  close  proximity  between  the  chief 
and  deputy,  equivalent  drag  effects,  and  the  only  perturbation  to  Keplerian  motion  is 
the  J2  effect. 

3.1.2  Yamanaka-Ankersen.  The  Yamanaka-Ankersen  (YA)  model  provides 
a  state  transition  matrix  mapping  Cartesian  initial  conditions  to  a  current  Cartesian 
state  [6] .  The  model  remains  as  a  hrst  order  linearization,  but  presents  a  closed  form 
solution.  Following  a  derivation  similar  to  the  HCW  construction,  the  differential 
equations  of  relative  motion  as  presented  in  the  YA  verbiage  are 

X  —koo^^'^x  +  2ujz  +  (jjz  +  u‘^x 

y  =  —ku^/'^y  +  Hf  +  acd  +  add  (3-8) 

^  2ku'^/‘^z  —  2ujx  —  (jjx  +  oj'^z 

where  uj  is  the  angular  rate  of  the  chief,  x,  y,  z  are  the  Cartesian  relative  motion 
coordinates,  k  =  aj  is  an  applied  force,  and  acd,Sidd  are  the  perturbing  forces 

on  the  chief  and  deputy.  The  Cartesian  coordinates  differ  from  the  reference  frame 
explained  in  Section  2.2.2.  The  previous  development  has  x,y,z  representing  the 
radial,  in-track,  and  cross-track  motion;  compared  to  this,  the  YA  model  places  its 
x,  y,  z  as  the  in-track,  negative  cross-track,  and  negative  radial.  Later,  matrices  of 
I’s  and  O’s  will  be  used  for  model  agreement.  Assuming  unforced  and  unperturbed 
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conditions,  the  differential  equations  can  be  simplified  as 


px"  —  2e  sin  9x'  —  e  cos  9x  =  2pz'  —  2e  sin  9z 

py"  —  2esi\i9y'  =  —y  (3.9) 

pz'  —  2e  sin  9z  —  (3  +  e  cos  9)z  =  —2px  +  2e  sin  9x 

where  the  prime  superscript  indicates  differentiation  with  respect  to  the  true  anomaly 
of  the  chief,  9  is  the  true  anomaly  of  the  chief,  and  p  =  1  +  e  cos  9.  A  solution  is  derived 
that  contains  the  same  integral  as  Eq.  2.2,  which  was  historically  left  in  integral  form 
until  the  Yamanaka-Ankersen  development.  The  YA  proposal  is  that  constant  angular 
momentum  allows  the  time  rate  of  change  of  the  true  anomaly  to  be  expressed  as 

(3.10) 

which  eventually  yields  the  expression 

k\t-to)  =  J{9)  (3.11) 

where  J{9)  is  the  integral  from  Eq.  2.2.  Following  a  geometric  interpretation  of  the 
relative  motion  and  some  matrix  algebra  to  prove  existence  and  uniqueness,  a  solu¬ 
tion  algorithm  is  proposed.  Beginning  with  initial  conditions  in  the  LVLH  frame,  a 
transformation  to  a  ’’tilde”  space  is  made.  The  tilde  space  is  non-physical  and  is  trans¬ 
formed  from  LVLH  space  through  multiplication  by  the  eccentricity,  true  anomaly, 
and  the  integral  J.  These  initial  tilde  conditions  are  then  transformed  into  a  set  of 
’’psuedoinitial”  conditions  that  can  then  be  propagated  forward  in  the  same  space. 
The  ’’psuedoinitial”  conditions  are  found  by  evaluating  the  inverse  of  the  STM  eval¬ 
uated  at  the  initial  time.  These  conditions  can  then  be  propagated  forward  in  tilde 
space.  The  out-of-plane  motion  does  not  require  the  ’’psuedoinitial”  transformation 
and  is  propagated  solely  in  the  tilde  state. 
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Assumptions  inherent  in  the  YA  model  are  close  proximity  between  the  chief 
and  deputy,  equivalent  drag  effects,  and  unperturbed  Keplerian  motion. 

3.1.3  Schweighart-Sedwick. 

3. 1.3.1  Model  Description.  The  Schweighart-Sedwick  (SS)  model  [36] 
is  perhaps  one  of  the  most  applicable  modes  of  this  study.  The  SS  system  assumes 
a  circular  reference  orbit,  just  as  the  assumption  inherent  in  the  ROE  derivation. 
However,  an  orbit  average  of  the  J2  perturbation  is  included  as  a  perturbing  force. 
The  model  is  presented  in  closed  form  with  initial  conditions  enforced  to  bound  the 
orbit.  Therefore,  the  SS  model  presents  a  method  of  initialization  for  a  periodic  orbit. 

The  derivation  of  the  equations  of  motion  is  rather  intuitive.  The  J2  effect  is 
inserted  as  an  acceleration  to  a  Newtonian  second  law  formulation.  Higher  order  terms 
are  linearized  with  respect  to  the  reference  orbit,  and  a  priori  corrections  are  made  to 
the  reference  orbit  to  account  for  the  perturbation  (including  period  matching,  nodal 
drift  correction,  and  cross-track  motion  correction).  The  model  uses  this  to  derive 
a  set  of  constant  coefficient,  linear,  ordinary  differential  equations  that  have  a  form 
rather  similar  to  the  HCW  differential  equations.  The  solutions  presented  in  [37]  are 
given  in  closed  form.  However,  the  desire  is  to  express  the  ROEs  as  a  function  of  any 
initial  condition.  The  differential  equations  are  derived  as 

X  —  2ncy  —  (5c^  —  2)n^x  =  0 

ij  +  2ncx  =  0  (3T2) 

^  -f  (3c^  —  2)n'^z  =  0 

where  x,  y,  z  are  the  radial,  in-track,  and  cross-track  coordinates,  n  is  the  mean  motion 
of  the  chief,  and  c  is  a  constant  incorporating  the  J2  effect.  For  this  study,  Eq.  3.12 
will  be  solved  in  closed  form  without  enforcing  initial  conditions.  The  SS  model 
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provides  a  periodic  orbit  of  the  following  form  (from  [37]) 


X  =  xq  cos  6  +  —yo  sin  9 
2c 

2c 

y  = - Xo  sm  9  +  yo  cos  9 

9 

z  =  sin(Q;  —  7) 

=  $0  cos  7o  sec  7 
7  =  tan“^(e  +  tan  70) 


(3.13) 


where  6*  is  a  modihed  time- varying  in-plane  phasing,  a  is  a  modihed  orbital  frequency, 
$  is  related  to  the  cross-track  amplihcation,  7  is  a  cross-track  phasing,  and  e  is  a 
frequency  related  to  the  change  in  cross-track  phasing.  The  various  constants  and 
coefficient  dehnitions  can  be  found  in  [37].  The  values  are  primarily  functions  of  the 
reference  orbit  inclination  and  are  first  order  in  J2.  The  in-plane  phasing  and  initial 
conditions  of  HCW  are  slightly  dampened,  and  the  cross-track  motion  operates  as  a 
function  of  angular  differences  rather  than  the  simple,  uncoupled  harmonic  oscillator. 
The  differential  equations  in  Eq.  3.12  will  be  solved  later. 

Assumptions  inherent  in  the  SS  model  are  close  proximity  between  the  chief  and 
deputy,  equivalent  drag  effects,  a  circular  chief  or  orbit,  and  the  only  perturbation  to 
Keplerian  dynamics  is  the  J2  effect. 


3.1.4  Numerically  Integrated  Truth  Model.  To  construct  a  foundation  for 
comparison  of  the  true  motion  for  both  Keplerian  and  perturbed  motion,  a  numerical 
integration  is  employed  in  Matlab  using  the  ode45  function  (a  one-step  solution  using 
the  Dormand-Prince  pair  of  the  Runga-Kutta  family  of  integrators).  Initial  inputs 
selected  are  the  chief  orbital  elements  [a,  e,  i,  u,  hi,  M]^  and  the  initial  LVLH  displace¬ 
ment  {p)  and  velocity  (p).  This  section  will  describe  the  equations  of  motion  for  the  J2 
case,  and  it  is  noted  that  the  unperturbed  motion  can  be  simply  obtained  by  setting 
J2  to  zero.  Given  the  chief  satellite  orbital  elements,  a  standard  two-body  formula¬ 
tion  will  yield  the  inertial  position  and  velocity  of  the  chief.  The  inertial  position  and 
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velocity  of  the  deputy  are  then  found  through  the  following  rotation 

Vrf  =  R*V+V-;  (3.14) 

where  the  rotation  matrix  rotates  a  vector  from  the  LVLH  or  orbital  frame  to  the 
inertial  frame,  and  is  dehned  as 

R  =  Or  00  Oh  (3.15) 

whose  column-vector  components  have  been  previously  dehned  in  Eq.  2.3.  The 
inertial  velocity  of  the  deputy  requires  a  Coriolis  term  in  the  rotation  and  is  found  as 

fd  =  h  +  R“  (/?+  X  p)  (3.16) 

with  /  having  been  previously  dehned  as  the  chief  true  anomaly,  and  co“  is  the  angular 
velocity  of  the  orbital  LVLH  frame  with  respect  to  the  inertial  frame.  Now  having  the 
inertial  position  and  velocity  of  the  chief  and  deputy,  initial  conditions  are  known  for 
the  integration.  Simultaneous  integration  of  Eq.  3.17  for  both  the  chief  and  deputy 
allow  for  their  states  to  be  known  at  any  time. 

(3.17) 

where  R(r')  is  the  J2  disturbing  function.  Knowing  the  state  at  a  time  t,  the  relative 
position  vector  can  be  found  as 

p  =  R“(f<,-r;)  (3.18) 

•  T 

where  R°*  =  R*°  and  rotates  from  the  inertial  to  the  orbital  LVLH  frame.  The 
relative  velocity  requires  a  cross  term  and  is  found  as 

p  =  (fd  -  fc  +  X  {fd  -  fc)^  (3.19) 
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The  J2  effect,  however,  modihes  the  angular  velocity  of  the  LVLH  frame  slightly. 
Rather  than  just  having  its  three-component  as  the  time  rate  of  the  true  anomaly, 
the  angular  velocity  using  mean  orbital  elements  is  found  from  Gim  and  Alfriend  [8] 
as 


/c 


A 


Q  sin  9  sin  i 
a)""  =  I  G  cos  6  sin  i 
9  +  Cl  cos  i 

where  the  mean  orbital  element  time  rates  are  detailed  in  Eq.  2.69. 


V' 


(3.20) 


7 


Note  that  as  expected  for  unperturbed  motion,  the  time  rates  of  the  perturbed 
orbital  elements  are  nulled  and  the  angular  velocity  vector  becomes 


00 


/o\ 

0 


and  Eq.  3.17  reduces  to 


and  the  previously  mentioned  rotations  remain  the  same. 


(3.21) 


(3.22) 


If  it  is  desired  to  examine  the  relative  motion  outside  the  linearized  regime 
the  analytical  models  allow,  the  numerically  integrated  truth  model  allows  for  the 
selection  of  initial  orbital  elements  for  both  the  chief  and  deputy  and  the  forward 
propagation. 


3.2  Methodology 

The  following  section  details  how  the  relative  motion  models  are  applied  to  the 
three  phases  of  this  study.  The  general  transformations  and  propagations  to  derive  the 
ROE  expressions  are  described.  Applications  of  ROEs  are  then  described  to  include 
initialization  of  a  bounded  relative  orbit  based  on  ROEs,  and  a  proposed  guidance 
and  navigation  algorithm. 
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3.2.1  Analytical  Expressions  for  the  Relative  Orbit  Elements.  Currently, 
expressions  for  the  ROEs  are  known  for  the  unperturbed  circular  chief.  It  is  now 
desired  to  examine  the  behavior  of  the  ROEs  for  a  perturbed,  circular  chief  and 
an  unperturbed,  elliptical  chief.  There  is  a  further  sub-division  in  this  study  between 
Cartesian  and  geometric  representations.  Let  the  transformation  from  Cartesian  space 
to  ROE  space  be  denoted  as 


R{t) 


ac(t)  x4t)  Vdit)  I3{t)  ^max  it)  %lj{t) 


=  A(X) 


(3.23) 


where  X  is  the  Cartesian  relative  position  and  velocity  states,  and  A  represents  a 
series  of  non-linear  operations  on  the  states.  The  analytical  instantaneous  states 
of  the  Yamanaka-Ankersen  and  Schweighart-Sedwick  models  will  be  substituted  as 
the  argument  in  the  A  transformation.  The  inverse  transformation  is  now  dehned 
as  X  =  A“^R(t).  Consequently,  this  allows  propagation  of  the  initial  relative  orbit 
elements  to  any  time  following  the  transformation.  The  direct  state  substitution  is 
employed  for  the  unperturbed  elliptical  chief  using  the  YA  model;  however,  the  cross¬ 
track  motion  of  the  Schweighart-Sedwick  motion  allows  for  a  more  intuitive  expression 
for  the  out-of-plane  relative  orbit  elements. 

Expressing  the  states  in  terms  of  orbital  element  differences  is  a  study  investi¬ 
gated  solely  for  the  perturbed  circular  chief.  The  primary  reasoning  for  this  restriction 
is  for  the  perturbed  or  two-body  elliptical  chief,  the  ROE  expressions  become  exceed¬ 
ingly  unwieldy.  However,  if  the  initial  circular  assumption  is  made  in  the  linearized 
mapping  between  the  Hill  frame  and  orbital  element  differences  [30],  changes  in  the  re¬ 
sulting  expressions  resulting  from  J2  can  be  obtained  rather  elegantly.  The  perturbed 
orbital  element  differences  will  be  taken  from  the  Gim-Alfriend  model.  A  chain  rule 
approach  will  also  be  taken  to  examine  the  time  rates  of  the  circular  ROEs  under  the 
J2  effect. 
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3.2.2  Applications  of  Osculating  Relative  Orbit  Elements.  The  following 
sections  briefly  detail  the  methodology  in  describing  the  applications  of  osculating 
ROEs. 


3.2.2. 1  Initialization  of  Relative  Orbit  Elements  for  Periodic  Trajecto¬ 
ries.  Several  methods  exist  (for  example,  Sabol  [3],  Schaub  [30],  and  Sengupta  [38]) 
to  set  the  initial  conditions  of  a  formation  to  induce  a  stationary  trajectory.  The  def¬ 
inition  used  in  this  study  for  a  stationary  orbit  is  that  of  a  relative  trajectory  that 
repeats  itself.  For  ideal  Keplerian  motion,  the  fundamental  requirement  for  a  static 
relative  trajectory  is  equal  periods  for  each  body  in  the  formation.  This  constraint 
equates  to  (5a  =  0  for  unperturbed  motion.  Although  the  constraint  remains  the  same, 
chief  eccentricity  complicates  the  expression  when  converted  to  Cartesian  coordinates. 
Considering  the  J2  perturbed  environment,  movement  of  the  perigee  and  nodal  re¬ 
gression  modify  the  periodicity  of  the  orbit.  Schaub  proposes  a  algorithm  [30]  for  J2 
invariance  that  results  in  selection  of  either  the  differential  inclination,  eccentricity,  or 
semi-major  axes  between  the  chief  and  deputy,  and  the  initial  differential  ascending 
node,  perigee,  and  mean  anomalies  calculated.  Moreover,  because  of  the  hrst-order 
approximation  of  the  J2  invariance  algorithm,  there  is  slight  long-term  drift. 

Initial  conditions  for  the  unperturbed  elliptical  chief  will  be  transformed  from 
the  Cartesian  LVLH  frame  to  ROE  conditions  using  an  equal  periodicity  condi¬ 
tion.  The  conditions  will  be  validated  using  the  Yamanaka-Ankersen  derived  ROEs. 
Bounding  the  perturbed  circular  orbit  will  use  the  Cartesian  conditions  set  in  the 
Schweighart-Sedwick  model  converted  to  ROE  space.  Validation  will  be  used  via  the 
Schweighart-Sedwick  derived  ROEs. 

3. 2. 2. 2  Guidance  and  Navigation.  Estimation  of  relative  orbit  ele¬ 
ments  also  carries  with  it  the  connotation  of  guidance.  As  a  means  of  showing  the 
applicability  of  the  ROEs,  the  eccentric  ROE  derivation  will  be  used  to  determine  an 
instantaneous  impulse  that  can  be  used  to  control  docking,  rendezvous,  and  guidance. 
Expressions  will  be  derived  to  numerically  solve  a  set  of  non-linear  equations  for  the 
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necessary  velocity  change  to  yield  a  desired  end  state.  Assuming  that  the  desired 
maneuver  is  a  solvable  problem,  the  calculated  post-burn  ROEs  will  be  compared 
with  the  desired  maneuver,  along  with  the  resulting  trajectories. 
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IV.  Osculating  Relative  Orbit  Elements 

The  following  sections  detail  the  mathematical  derivations  necessary  to  obtain  analyt¬ 
ical  expressions  for  the  relative  orbit  elements  with  an  eccentric  chief  and  a  perturbed 
circular  chief.  Section  4.1  details  the  derivation  for  the  elliptical  chief.  Section  4.2  de¬ 
rives  the  closed  form  expressions  of  the  J2  perturbed  osculating  ROEs.  Finally,  Sec. 4. 3 
details  areas  of  applications  for  the  ROEs  using  a  guidance  and  navigation  algorithm 
and,  in  addition,  analytical  expressions  to  initialize  a  bounded  relative  orbit. 


4.1  The  Non- Circular  Chief 

4.1.1  State  Propagation.  Yamanaka  and  Ankersen  presents  an  algorithm  for 
calculating  relative  satellite  motion  [6] .  Before  detailing  the  matrix  multiplication,  the 
verbiage  used  in  the  YA  model  is  now  introduced  as 


p  =  1  +  e  cos  6 
s  =  psinO 
c  =  p  cos  9 
s'  =  cos  9  +  e  cos  29 
S  =  —(sin  6^  -I-  e  sin  29) 


J 


h 


pZ 


{t  -  to) 


(4.1) 


Here,  9  is  the  true  anomaly  of  the  chief,  p  is  the  semi-latus  rectum  of  the  chief’s  orbit, 
and  h  is  the  chief  orbital  angular  momentum.  Given  in  [6],  the  in-plane  motion  can 
be  propagated  as 


X 

1  (1  +  p)  -s  (1  +  j) 

Xo 

z 

Os  c  (2  —  3esJ) 

zo 

Vx 

0  2s  2c  —  e  3(1  —  2esJ) 

^xO 

Vz 

0  s'  c'  — 3e(s'J  -I- 

VzO 

(4.2) 
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It  is  very  important  to  note  that  the  terms  x^z^v^^Vy  terms  are  reported  here  using 
the  YA  LVLH  frame.  A  conversion  will  be  made  later  for  model  agreement.  For 
the  moment,  the  4x4  matrix  in  Eq.  4.2  is  denoted  as  <F,  such  that  Xip  =  ^XipQ, 
with  the  in-plane  YA  states  used  as  arguments  in  the  unbolded  vector  X^p.  The  tilde 
and  bar  accents  denote  transformed  coordinates;  however,  the  transformation  is  not 
a  rotation  but  can  be  shown  to  be  a  product  of  matrix  multiplication.  Letting  the 
Cartesian  to  tilde  space  transformation  be  denoted  by  H,  and  the  tilde  to  Cartesian 
transformation  as  H~^,  the  matrices  are  found  by  inspection  of  the  YA  model  as 


such  that 


H 


H-i 


P  0 

0  p 

—es  0 


0  0 
0  0 


0 


—es  0 


1 

k'^p 


i  0  0  0 

p 

0^00 

p 

k'^es  0  k‘^p  0 


0  k'^es  0  k'^p 


X  =  HX 
X  =  U^^X 


(4.3) 


(4.4) 


The  barred  coordinates  are  a  function  of  the  tilde  coordinates.  Allowing  the  trans¬ 
formation  denoted  by  L  such  that  X  =  LX,  the  matrix  L  is  expressed  as 


L 


l-e^ 


1 

1  -e2 


0 

0 


0 


3p  -I-  —  1 


—es  ^  j  — ec  -|-  2 
s{l  +  l)  c-2e 

c(l  +  i)+e  -s 

—p^  es 


(4.5) 
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The  final  matrix  multiplication  for  the  in-plane  motion  then  becomes 


X,p{t)  =  li-\t)X,p{t) 

Xip{t)  =  U~^{t)^{t,to)Xipo 

Xip{t)  =  U~^{t)^{t,to)L{to)Xipo 
X,p{t)  =  H-\t)<l>{t,to)L{to)H{to)X,po 


(4.6) 


There  now  exists  a  single  state  transition  matrix  to  propagate  the  initial  conditions 
to  a  desired  time  for  the  in-plane  motion.  Using  the  YA  defined  STM,  the  out  of 
plane  motion  is  found  more  simply  as  a  propagation  in  YA  tilde  coordinates.  Letting 
the  out-of-plane  tilde  to  Cartesian  transformation  be  denoted  as  Hi,  the  propagation 
goes  as 

ATop  =  ^(t)<hop(f,  to)Hi(to)Yopo  (4.7) 

with  the  matrices  defined  as 


^op(t,to) 

Hr'(f) 

Hi  (to) 


c  s 
—s  c 


Pe-do 


e-do 

1 

p 


0 


k‘^es  k'^p 


P  0 
1 


-es 


k'^pA  e- 


e-eo 


(4.8) 


where  9—6q  indicates  that  values  in  the  matrix  are  calculated  at  the  difference  between 
the  initial  and  current  true  anomaly  of  the  chief.  However,  the  Yamanaka-Ankersen 
model  uses  a  coordinate  system  that  does  not  agree  with  the  nomenclature  already 
established  in  the  current  study.  The  YA  model  labels  its  x,  y,  z  coordinates  as  in¬ 
track,  negative  cross-track,  and  negative  radial,  counter  to  the  current  study’s  usage 
of  x,y,z  as  radial,  in-track,  cross-track  (Xjioe)-  This  requires  a  multiplication  of 
the  initial  conditions  and  end-state  by  a  matrix  of  ones  and  zeros.  Labeling  this 
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’agreement’  matrix  as  U  such  that  X^qe  =  ,  the  in-plane  and  out-of-plane 

portions  are  found  by  inspection  as 


0-100 
0  0  0  -1 
10  0  0 
0  0  10 

-1  0 
0  -1 


Uop 


This  conversion  in  matrix  form  is  now  seen  as 


X 

0 

-1  0 

0 

y 

0 

0  0 

X 

1 

0  0 

0 

y 

ROE 

0 

0 

0 

z 

-1 

0 

y 

z 

ROE 

0 

-1 

y 

VA 


(4.9) 


(4.10) 


The  hnal  propagation  in  the  nomenclature  of  the  current  study  now  goes  as 


=  U,pH-\tMt,to)L(to)H(toM^X,,(to) 

Xop(t)  =  UopH{\t)^op(t,to)Hi(to)U^pXopo 


(4.11) 


where  Xip(t)  = 


iT 


z  z 


.  The  X,  y,  z  triplet  now  denotes 


X  y  X  and  Xop{t)  — 

radial,  in-track,  and  cross-track  coordinates.  Condensing  Eq.  4.11,  the  expression 
can  be  written  as 

Xip{t)  =  ^ip{t,to)Xip{to) 

Xopit)  —  ^  ip  (^7  to)^OpO 


(4.12) 
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where  to)  is  a  4x4  matrix  propagating  the  in-plane  trajectory,  and  4>op(t,  to)  is 
a  2x2  matrix  propagating  the  out-of-plane  trajectory.  Now  having  the  state  propa¬ 
gation,  the  relative  orbit  elements  can  now  be  solved  for  analytically. 


4- 1-2  Analytical  Derivation  of  Relative  Orbit  Elements  for  the  Unperturbed 
Noncircular  Chief.  Now  we  can  derive  the  ROEs  for  the  noncircular  chief  using 
the  full  matrix  multiplication  of  Eq.  4.11.  For  now,  the  derivation  will  not  make  any 
assumptions  on  initial  conditions  (In  Section  4.1.3,  a  perigee  epoch  assumption  will 
be  enforced  and  justihed).  To  avoid  cumbersome  equations,  the  components  of  the 
state  transition  matrices  in  Eq.  4.11  will  be  expressed  for  the  in-plane  motion  as 


X 

All  •  •  • 

.  .  .  Ai4 

Xq 

y 

yo 

X 

Xo 

y 

A41  . . . 

. . .  A44 

yo 

and  for  the  out  of  plane  motion  as 


z 

Cu 

C12 

Zo 

z 

C21 

C22 

Zo 

(4.13) 


(4.14) 


The  analytical  expressions  for  the  parameters  Rii,Ri2, . . .  ,6*22  can  be  found  in  Ap¬ 
pendix  C.  However,  it  is  important  to  note  that  the  full  expressions  are  very  cumber¬ 
some  and  are  extremely  tedious  to  present  analytically.  Those  components  presented 
in  the  Appendix  utilize  an  assumption  of  a  perigee  initial  condition  and  epoch  time 
of  to  =  0.  This  is  not  without  loss  of  applicability  as  is  noted  in  Section  4.1.3.  The 
following  derivations  are  performed  using  generic  Aij  components  and  is  valid  for  both 
the  full  and  assumed  expressions. 
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4- 1.2.1  Initial  ROE  Propagation  in  lieu  of  the  State  Vector. 
in-plane  initial  conditions  expressed  in  terms  of  ROEs  are 


-^cos/3o  +  a;do 

Oeo  sin/3o  +  ym 

sin  /3o 

Oeon  cos  Po  -  \nXdo 


(4.15) 


With  respect  to  the  ROEs,  the  non-phasing  terms  can  be  separated  from  the  angular 
terms  by  inspection  as 


-^cos 


+  XdO 


30  sinpQ  +  ydo 

^nsiii/do 


i  cos  Po  ■ 


sin  /3o  0  1 

^  sin  /3o  0  0 


Inxdo 


n  cos  Po  — 


—  0 
2  ^ 


(4.16) 


Rewriting  the  initial  Cartesian  vector  as  Eq.  4.16  and  substituting  in  Eq.  4.13,  the 
in-plane  state  can  be  written  as 

xl  r (Ri2  sin  Po  -  +  Rl4n  cos  Po  +  ^ 

cigo 

y  _  {A22smpo  -  +  A24ncospo  +  [A21  -  ^)  R22 

X  (R32Sin/3o  -  +  AsiUCOSpo  +  (^ASl-^)  R32 

y\  [(R42sin/3o  -  +  A,^ncosPo  +  ^aAI-^)  R42J 

(4.17) 


For  simplicity,  each  component  of  Eq.  4 
the  time  varying  matrix  B(t)  such  that 


.  4.17  will  be  assigned  to  parameter  Bij  to  form 


=  B(t)  XdO 


(4.18) 
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Similar  to  the  initial  set-up  of  the  in-plane  motion,  the  out-of-plane  state  can 


be  expressed  as  a  function  of  the  initial  ROEs  using  the  conversion 


Zo  = 

Zo  =  ZmaxOn  COS 


(4.19) 


Substituting  Eq.  4.19  in  Eq.  4.14,  the  cross-track  state  can  be  found  as 


z 
z 

Factoring  out  the  Zmaxo  term  and  multiplying  the  resulting  inner  matrices,  a  single 
multiplication  can  be  found  as 


On 

O12 

C21 

1 

(N 

z^axo  sin^/>o 

ZmaxO'l^  COS  i/Jq 


(4.20) 


Cii  sin  ipQ  +  Ci2n  cos  ipQ 
C21  sin  il)Q  +  C22n  cos  ^jJo 


^maxO 


(4.21) 


and  assigning  each  component  of  the  matrix  in  Eq.  4.21  the  parameter  Dij,  the 
cross-track  motion  is  now  propagated  as 


2: 

i 


D(^)^maa:0 


(4.22) 


4- 1-2.2  State  Substitution  as  Arguments  in  the  Relative  Orbit  Element 
Expressions.  Having  expressions  for  the  state  as  a  function  of  initial  ROEs,  the 
propagated  ROEs  can  now  be  found.  The  relative  semi-major  axis  is  found  as 


Cte 


(4.23) 


Expanding  the  above  expression  becomes  a  rather  tedious  process.  To  gather  enough 
information  to  observe  the  impact  of  the  coefficients  on  the  initial  conditions,  the 
expression  for  Oe  is  squared,  expanded,  and  grouped  according  to  the  initial  ROEs. 
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Foregoing  the  extensive  algebra,  the  highly  coupled  end  result  is  found  as 


o-e  —  y  C^l^eO  +  ^'2^‘dO  +  ^^VdO  +  CX^^doVdO  +  +  Ct6®e0l/d0  (4.24) 

where  the  values  for  the  time-varying  a  coefficients  are  listed  in  Appendix  C.  More¬ 
over,  the  derivations  are  given  in  further  detail  in  Appendix  B. 

In  contrast  to  the  relative  semi-major  axis,  the  radial  displacement  is  a  linear 
combination  of  two  Cartesian  states.  Expressing  the  radial  displacement  as 

Xd  =  4:X  +  2-  (4.25) 

n 

state  substitution  yields  after  grouping  like  terms 


Xd  —  o'laeo  +  o'2Xdo  +  cr^ydo 


(4.26) 


where  the  values  for  the  time- varying  a  coefficients  are  listed  in  Appendix  C. 

Very  similar  to  the  radial-displacement  and  expressing  yd  SiS  i/d  =  y  —  ^  ,  the 
in-track  displacement  is  also  found  as  a  function  of  time-varying  coefficients  as 


Ud  —  SiOeO  +  ^2XdO  +  ^SUdO 


(4.27) 


with  the  S  expressions  listed  in  Appendix  C. 

The  in-plane  phasing  angle  is  expressed  as 


f3  =  tan 


-1 


X 


3nx  -|-  2y 


(4.28) 


which  in  reality  is  the  ratio  of  the  two  parenthetical  components  of  the  semi-major 
axis  function.  The  in-plane  phasing  is  found  as 


f3  =  tan 


-1 


-BsifleO  +  B32XdO  +  B^^ydO 


{3nBii  +  2i?4i)  Oeo  +  {3nBi2  +  2B42)  Xdo  +  (3ni?i3  -|-  2B43)  ydo 


(4.29) 
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where  the  B  coefficients  have  previously  been  dehned  from  Eq.  4.17. 


The  out-of-plane  phasing  angle  ifj  is  nonlinear  in  the  arctangent  function.  The 
expression  for  ip  is 

Ip  =  tan“^  (4.30) 

and  after  substitution  it  is  found  as 


where  the  Dij  terms  are  obtained  from  Eq.  4.22.  The  hnal  ROE  to  examine  is  the 
cross-track  amplihcation.  Using  the  conversion 

(1)0.-  (4.32) 

which  is  found  as  a  time- varying  value  as 


D 


21 


+  /1? 


11 


^max  ^maxO 


(4.33) 


To  provide  a  summarization  of  the  osculating  ROEs  with  an  eccentric  chief,  the 


following  list  is  provided 


fle  —  a/  dlOgQ  +  0'2x'^q  +  ttal/rfo  (^i^dOUdO  +  +  tteOeOl/dO 


p  1 

r 

1 

Q-eO 

Xd 

= 

0-1 

<72 

<73 

XdO 

Vd 

Si 

S2 

S3_ 

UdO 

f3  =  tan  ^ 
Ip  =  tan 


Bsitteo  +  B^iXdo  +  B^^ydo 


{SnBii  +  2R41)  Oeo  +  (3n-Bi2  +  2R42)  +  (SnRis  +  2R43) 

/ nDii 


21 


^max  XynaxO  ^ 


D 


21 


n 


+  D^ 


11 


(4.34) 


with  the  time- varying  coefficients  q;*,  Ui,  Sj,  B^j,  and  Dij  detailed  in  Appendix  C. 


4.1.3  Perigee  Epoch  Simplification.  The  time- varying  coefficients  of  Eq. 
4.34  are  extremely  cumbersome.  To  simplify  the  complicated  coefficients,  it  is  now 
assumed  that  the  epoch  time  is  equivalently  zero  (to  =  0),  and  that  the  true  anomaly 
of  the  chief  at  epoch  is  perigee  (6^0  =  0).  The  implication  of  this  assumption  is  not  in 
immediate  disfavor  of  the  model.  From  a  deterministic  standpoint  and  employment 
of  a  state  transition  matrix,  knowledge  of  the  state  at  any  one  time  allows  knowledge 
of  the  state  any  other  time.  For  example,  consider  a  set  of  relative  navigation  data  at 
a  time  t  that  yields  a  relative  state  X(t).  From  the  fundamentals  of  state  transition 
matrices  given  in  Appendix  A,  the  state  at  perigee  at  time  tp  can  be  found  as  X  (tp)  = 
^fip,t)X(t)  =  ^~^(t,tp)X(t),  and  the  model  reinitialized  to  t  =  0  at  the  perigee 
time  for  model  agreement  presents  little  difficulty.  The  only  spatial  limitations  on 
the  initial  conditions  are  those  limited  by  the  degree  of  accuracy  in  the  linearization 
assumptions  in  the  model. 

The  actual  enforcement  of  the  initial  condition  assumption  allows  for  the  L  and 
H  of  Eq.  4.4  and  Eq.  4.5  transformations  to  greatly  simplify.  However,  the  nature 


of  the  derivation  demands  that  the  general  form  of  Eq.  4.34  remains  the  same.  The 
coefficients  simplify  signihcantly  and  are  detailed  in  Appendix  C. 

4- 1-4  Numerical  Examples.  To  demonstrate  the  applicability  of  the  oscu¬ 
lating  parameterization,  the  following  numerical  examples  will  demonstrate  the  accu¬ 
racy  of  the  method  when  applied  to  numerical  integration  of  two-body  (Keplerian) 
motion.  The  relative  states  are  calculated  analytically  using  the  ROE  parameteriza¬ 
tion  and  compared  to  the  resulting  integration.  Stationary  (repeating  relative  orbits) 
examples  of  varying  eccentricity  are  given,  followed  by  a  drifting  example.  The  sta¬ 
tionary  relative  orbits  are  found  using  the  conditions  detailed  in  Appendix  F.  Initial 
conditions  are  given  in  the  form  Rq  =  Oeo  Xdo  Vdo  Po  Zmaxo  i’o  units  of 

(km,rad)  where  appropriate,  and  the  chief  orbital  elements  are  given  in  the  form 
Gco  =  a  e  i  u  Mq  with  units  of  (km,rad)  as  well. 

4-1  ■4-1  Low-Eccentricity,  Stationary  Relative  Orbit.  Given  the  follow¬ 
ing  initial  conditions 

Ro  =  [3.4890  0.1  0.5  0  0.5  0 

^  ,  (4.35) 

eco=  10000  0.01  0.5236  0  0  0 

The  three-dimensional  propagation  using  osculating  ROEs  as  compared  to  the  integra¬ 
tion  is  shown  in  Fig.  4.1  By  immediate  inspection,  the  three-dimensional  trajectory 
using  the  YA  dynamics  compares  well  with  the  numerical  integration.  Little  skewness 
is  observed  in  the  out-of-plane  motion,  and  the  projected  in-plane  resembles  a  closed 
irregular  conic. 

Figure  4.2  provides  a  display  of  the  in-plane  ROEs  as  calculated  through  the 
derived  equations  compared  to  the  HCW  prediction.  The  x-axis  is  labeled  as  orbit 
fraction  and  represents  the  current  time  in  the  orbit  divided  by  the  nominal  chief 
period.  The  ROEs  for  a  circular  chief  predict  a  constant  Og  and  Xd]  however,  using 
the  eccentric  chief  parameterization  shows  a  periodic  Og  that  coincidentally  intersects 
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Relative  Trajectories 


YA 

2BP 


Figure  4.1:  Three-Dimensional  Relative  Trajectory  using  Osculating  ROEs  Compared 
to  Numerical  Integration  with  Chief  Eccentricity  of  0.01 

the  HCW  semi-major  axis  at  perigee.  The  same  condition  is  true  for  the  radial 
displacement  with  a  HCW  coexistence  at  perigee.  Using  the  HCW  realization  would 
introduce  drift  in  this  case,  but  the  eccentricity  effects  yield  a  bounded  trajectory. 
An  energy  argument  logically  provides  that  the  periods  are  no  longer  matched  with 
the  added  eccentricity.  The  in-plane  phasing  shows  little  deviation  from  the  linear 
prediction  of  HCW. 

Figure  4.3  provides  the  out-of-plane  ROEs;  the  amplitude  is  slightly  perturbed, 
reaching  a  minimum  at  apogee,  and  osculating  the  HCW  value  at  perigee.  The 
maximum  out-of-plane  error  between  the  HCW  realization  and  the  eccentric  chief 
value  is  on  the  order  of  0.01  km.  These  observations  imply  that  using  the  circular 
ROEs  for  low  eccentricity  may  be  a  favorable  trade-off  between  model  complexity  and 
accuracy. 
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Orbit  Fraction 


Figure  4.2:  In-Plane  Osculating  Relative  Orbit  Elements  with  Chief  Eccentricity  of 
0.01 


Figure  4.3:  Out-of-Plane  Osculating  Relative  Orbit  Elements  with  Chief  Eccentricity 
of  0.01 
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4.1. 4-2  Medium- Eccentricity,  Stationary  Relative  Orbit. 


Given  the 


following  initial  conditions 


Ro  — 


1.3462  0.5  0.5  0  0.5  0 


GcO  — 


12000  0.3  0.5236  0  0  0 


(4.36) 


The  three-dimensional  propagation  using  osculating  ROEs  as  compared  to  the  inte¬ 
gration  is  shown  in  Fig.  4.4  The  propagation  compared  with  the  numerical  integration 

Relative  Trajectories 

- YA 

- 2BP 


Figure  4.4:  Three-Dimensional  Relative  Trajectory  using  Osculating  ROEs  Compared 
to  Numerical  Integration  with  Chief  Eccentricity  of  0.3 


is  near  exact.  There  are  additional  frequencies  manifesting  in  the  cross-track  motion, 
skewing  the  trajectory.  The  projection  still  remains  a  closed  conic,  with  a  shape  that 
is  also  skewed  by  additional  harmonics. 

Figure  4.5  provides  the  in-plane  ROEs  compared  to  the  expected  HCW  values. 
Immediate  equality  is  noted  again  at  perigee  due  to  the  absence  of  radial  velocity. 
The  addition  of  eccentricity  accelerated  the  HCW  drift  rate  when  compared  to  the 
eccentric  motion.  The  semi-major  axis  now  becomes  more  of  an  indication  of  the 
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magnitude  of  the  radial  and  in-track  oscillations,  which  have  significant  less  phase 
than  the  circular  counterpart. 


Figure  4.5:  In-Plane  Osculating  Relative  Orbit  Elements  with  Chief  Eccentricity  of 
0.3 

The  eccentricity  effects  seen  in  Fig.  4.6  decrease  the  magnitude  of  the  cross¬ 
track  motion  on  the  order  of  0.25  km.  The  cross-track  phasing  lags  with  respect  to 
the  HCW,  showing  the  variance  in  the  orbital  angular  momentum  in  comparison  to 
the  circular  assumption. 
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Figure  4.6:  Out-of-Plane  Osculating  Relative  Orbit  Elements  with  Chief  Eccentricity 
of  0.3 


4- .1.5  Physical  Interpretation  of  the  Osculating  Relative  Orbit  Elements  for 
the  Unperturbed  Non-Circular  Chief.  For  the  circular  chief,  the  ROEs  elegantly 
describe  the  geometry  of  the  relative  motion.  When  applied  to  the  current  scenario, 
the  physical  signihcance  of  the  ROEs  encroaches  a  mathematical  abstract.  The  effects 
of  chief  eccentricity  were  detailed  in  Section  2.5.1,  but  can  also  be  observed  in  time- 
response  of  the  newly  derived  ROEs. 

Focusing  solely  on  the  radial  and  in-track  motion,  the  in  plane  projection  of  the 
relative  trajectory  is  found  now  as  the  locus  of  points  {x,  y)  such  that 


/ a:(t)  ^  / y{t)-yd{t)\  ^  ^ 

ae(f)  J  2a,{t)  J 


(4.37) 


This  implies  directly  that  the  position  {x,y)  is  the  solution  to  Eq.  4.37  given  the 
ROE  parameterization  R{t).  This  describes  the  motion  of  a  2x1  ellipse  whose  center 
is  translating  in  the  plane  according  to  the  trajectory  of  {xd,  yd)-  The  semi-major  axis 
of  the  instantaneous  ellipse  is  time- varying  according  to  the  function  ae{t).  Moreover, 
as  the  ellipse  may  intersect  with  the  true  trajectory  at  multiple  points,  the  location 
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described  at  time  t  can  be  located  through  the  phasing  angle  f3.  Figure  4.7  provides 
a  visualization  of  this  application,  where  E  denotes  the  osculating  ellipse,  and  R{t)  is 
the  ROE  parameterization  at  time  t.  This  description  hts  perfectly  with  the  idea  of 


Figure  4.7:  In-Plane  Relative  Orbit  Element  Characterization  for  the  Eccentric  Chief 

an  osculating  relative  orbit.  If  at  a  time  t,  the  calculated  ROE  space  were  initialized, 
the  osculating  ellipse  entirely  describes  the  relative  trajectory  if  the  chief  eccentricity 
were  not  present.  However,  the  inclusion  of  the  eccentricity  terms  yields  a  different 
response  on  the  true  trajectory,  and  the  osculating  ellipse  simply  coincides  with  the 
instantaneous  position  terms  as  a  function  of  the  parameterization.  An  interesting 
note  is  that  when  the  orbit  is  bounded,  as  the  chief  approaches  the  initial  conditions, 
the  osculating  ellipse  becomes  coplanar  with  the  HCW  ellipse. 

The  motion  of  the  osculating  ellipse,  termed  by  the  author  as  the  osculational 
translation,  is  also  an  indication  of  the  offset  of  the  ROE  parameterization  from  the 
true  solution.  The  osculating  semi-major  axis  scales  the  magnitude  and  places  the 
in-plane  components  with  respect  to  the  phasing,  and  the  offsets  compensate  for  the 
remainder  of  any  deviation.  A  typical  plot  of  the  osculational  translation  is  found  in 


65 


Fig.  4.8  where  the  (x^,  Vd)  trajectory  the  center  of  the  given  osculating  ellipse.  (This 


Figure  4.8:  Osculational  Translation 

example  was  generated  using  a  chief  with  9000  km  semi-major  axis,  0.3  eccentricity, 
15  degree  inclination,  and  all  other  orbital  elements  as  zero,  with  a  Cartesian  LVLH 
initial  condition  of  [p,  =  [0,  0.005,  0.001,  0.001,  0,  0.001]^(km,km/s)). 

The  cross-track  motion  is  far  less  complicated  for  this  unperturbed  case.  The 
out-of-plane  motion  remains  uncoupled  but  is  attenuated  from  the  HCW  solution  as 
a  function  of  the  eccentricity  as 

z{t)  =  ZmaxoZ{t)  (4.38) 

The  form  of  Eq.  4.38  shows  the  HCW  cross-track  maximum,  but  is  modihed  by  the 
product  Z{t)  sin 'll)  [t).  This  will  be  referred  to  as  the  cross-track  attenuation,  and 
forces  a  dampening  in  amplitude  and  skewness.  Figure  4.9  provides  an  example  of 
the  cross-track  motion  normalized  by  the  HCW  expected  maximum  (using  the  same 
initial  conditions  as  the  previous  example).  Inspection  of  Fig.  4.9  reveals  extrema 
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0.8, 


Figure  4.9:  Cross- Track  Attenuation 

values  of  the  attenuation,  showing  an  actual  maximum  amplihcation.  This  will  occur 
when  the  cross-track  phasing  is  orthogonal  to  the  in-plane  motion  {'ip  =  0).  This 
shows  a  highly-coupled  relationship  with  the  true-anomaly  and  mean  motion  of  the 
chief,  and  also  indicates  a  maxima  on  an  elliptical  cylinder  containing  the  entirety 
of  the  stationary  trajectory,  or  of  an  instantaneous  elliptic  cylinder  containing  the 
instantaneous  motion  of  the  drifting  trajectory.  In  addition,  the  cross-track  phasing 
simply  provides  directionality  to  the  magnitude. 

At  any  point  in  the  trajectory,  the  instantaneous  ellipse  shows  the  behavior  of 
the  HCW  realization.  This  means  that  if  the  model  were  initialized  at  any  time  t, 
the  osculating  ellipse  is  the  HCW  realization  using  the  relative  state  at  t  as  its  initial 
conditions.  It  has  been  demonstrated  that  the  eccentric  ROEs  equate  to  the  circular 
ROEs  at  perigee.  This  is  resulting  from  the  fact  that  for  any  closed  orbit,  the  velocity 
at  perigee  is  purely  tangential  with  no  radial  component.  At  each  other  point  in  the 
relative  trajectory,  the  chief  has  a  radial  component  that  skews  the  ideal  2x1  ellipse. 
The  idea  of  the  osculating  ellipse  becomes  a  mathematical  abstract  to  describe  the 
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resulting  relative  trajectory,  as  in  this  sense  the  geometric  parameterization  is  now 
a  locus  of  solutions  to  a  time  variant  ellipse  .  Although  the  physical  signihcance  of 
the  ROEs  tends  to  break  down,  the  parameters  still  describe  certain  properties.  The 
following  details  the  physical  properties  of  the  eccentric  ROEs 

•  The  parameter  Og  does  describe  the  semi-major  axis  of  the  in-plane  projection 
of  the  HCW  trajectory  if  the  model  were  initiated  at  at  time  t  with  initial  con¬ 
ditions  x{t),y{t).  However,  it  also  helps  determine  the  instantaneous  maximum 
amplitude  of  the  radial  (ae(t)/2)  and  in-track  motion  (ae(t)). 

•  The  parameters  Xd  and  yd  do  detail  the  center  of  the  osculating  ellipse.  However, 
in  the  eccentric  ROE  scenario,  they  describe  the  offset  of  the  radial  and  in-track 
motion. 

•  The  angles  (3  and  'ijj  remain  as  indications  of  the  in-plane  and  cross-track  phasing 
frequencies,  but  lose  true  physical  interpretation. 

•  The  parameter  Zmax  indicates  the  instantaneous  maximum  magnitude  of  the 
cross-track  motion. 


4-2  The  Perturbed,  Circular  Chief 

4.2.1  Relative  Orbit  Elements  for  the  Perturbed,  Cireular  Chief  Using  Geo- 
metrieal  Insight  and  Linearized  Mapping.  It  has  been  shown  that  relaxation  of  the 
circular  assumption  complicates  the  physical  interpretation  of  the  ROEs.  However, 
their  applicability  to  low  eccentricity  references  orbits  still  holds.  The  focus  of  this 
study  now  moves  to  how  the  equatorial  bulge  of  the  psuedo-spherical  Earth  impacts 
the  circular  development. 

4. 2. 1.1  Derivation  of  Relative  Orbit  Elements  for  the  Unperturbed  Cir¬ 
eular  Chief  using  Orbital  Element  Differenees.  As  a  means  to  gain  expectations  for 
the  behavior  of  the  relative  geometry,  an  analytical  exercise  is  hrst  performed,  using 
the  Gim-Alfriend  mapping  between  Cartesian  and  orbital  element  differences  [8] 

X(t)  =  {Ait)  +  aB{t))  6e{t)  (4.39) 


Where  Se{t)  =  —  Oc  are  the  orbital  element  differences  at  a  time  t.  To  gather  a 

more  intuitive  understanding,  the  restriction  is  made  to  work  in  mean  orbital  element 
space.  The  unperturbed  portion  A  from  condenses  to 


A  = 


10  0  —RcosO 

0  0  0  Id  sin  6* 

0  R  0  0 

0  0  2VtCos9 

0  0  RsinO  0 

0  0  Id  cos  9  0 


—Rsm9  0 

—Id  cos  6^  0 

0  R  cos  i 

2Vtsm9  0 

0  —Rsmicos9 
0  Id  sin  9 


(4.40) 
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while  the  mean  perturbed  B  portion  from  [8]  reduces  to 


B  = 


0 

0 

0 


7n  cos^  i 
0 


7n  cos  B  sin  i  cos  i 


0 

0 

0 

0 

0 

0 


0  0 

(^5  cos^  z— 1^— n  sin  0) 
^  4i? 

0  0 

n  sin  i  cos  i  Ci 

2R  ^ 

0  0 

2n  cos  0  sin^  i  o 

4i7  ^ 


0 

(^5  cos^ 

IR 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


(4.41) 


At  this  point,  it  is  desired  to  express  the  unperturbed  relative  orbit  elements  as  func¬ 
tions  of  the  unperturbed  orbital  element  differences;  this  is  achieved  by  setting  a  equal 
to  zero  in  Eq.  4.39.  Also  inserting  the  circular  assumption,  the  state  transformation 
becomes 


X 

1 

0 

0 

—a  cos  9 

— asin^ 

0 

5a 

X 

0 

0 

0 

na  sin  9 

—na  cos  9 

0 
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y 
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a 
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0 

0 

acosi 

5i 

y 

3n 

2 

0 

0 

2na  cos  9 
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0 

5qi 

z 

0 

0 

asin6' 

0 

0 

— a  sin  7  cos  6^ 

5q2 

z 

0 

0 

na  cos  9 

0 

0 

na  sin  9 

5Q 

(4.42) 


where  6  is  now  the  argument  of  true  latitude,  Vt  is  the  tangential  velocity  (mean 
motion  for  this  specific  example),  and  all  other  variables  are  previously  defined.  This 
can  be  expanded  and  substituted  into  the  ROE  expressions  with  easy  verification  to 
find 

Oe  =  2a\J  {Sqif  +  {Sq2f  =  2aed 
Xd  =  6a 


ijd  =  a56  +  aSQ  cos  i  —  2a  {5qi  sin  9  —  cos  65q2) 
{ 5qism.9  -  5q2Cos9' 

jj  —  tciri 


6qi  cos  9  +  Sq2  sin  9 

=  a\J  {5i)^  +  ((50sinf)^ 

5i  sin  9  —  SQ  cos  9  sin  i 


Ip  =  tan 


-1 


5i  cos  6^  -|-  (50  sin  9  sin  i 


(4.43) 
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Immediately,  it  is  noted  that  the  boundedness  condition  of  =  0  easily  extends  to 
orbital  element  different  space  as  (5a  =  0  (the  equal  period  condition).  Also  of  note 
is  that  the  non-constant  values  {yd,j3,'ip)  remain  the  same.  As  a  direct  fallout,  these 
expressions  are  near  identical  to  parameters  developed  by  Schaub  [30].  When  the 
expression  given  in  Eq.  4.43  for  Zmax  is  divided  by  the  semi-major  axis,  the  resulting 
expression  is  given  by  Schaub  as  the  angle  between  the  orbit  planes  of  the  deputy  and 
the  chief.  This  implies  that  the  Zmaxio.  is  a  spherical  tilt  angle  between  the  two  orbit 
planes  for  very  small  orbital  element  differences.  The  inquisitive  response  to  this  is 
how  the  J2  perturbation  will  affect  these  expressions. 

4-2. 1.2  Osculating  Relative  Orbit  Elements  from  the  J2  Effect  using  Or¬ 
bital  Element  Differences:  An  Analytical  Approach.  The  objective  now  is  to  derive 
the  mean  effects  of  J2  on  the  relative  orbit  elements.  Using  the  mean  secular  J2  time 
rates,  the  following  linear  relationships  are  developed  based  on  the  secular  drift  in  the 
orbital  elements 

5a  =  5ao 

<^(li  =  <li  =  (Ro  cos  (OAt)  —  q2o  sin  (OAt) 

=  Cd  cosu(t) 

5q2  =  q2  =  Qio  sin  {uAf)  +  ^20  cos  [uAf) 

(4.44) 

=  edAnuiif) 

5Vt  =  (5f2o  +  StlAt 
66  =  (56*0  -f  66 At 
6i  =  (5*0 

where  is  the  eccentricity  of  the  deputy.  Direct  state  substitution  can  now  be  used  to 
examine  the  effect  on  the  unperturbed  ROEs.  Beginning  with  the  relative  semi-major 
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axis 


tte  =  2a\J  {6qif  +  {6q2f 

=  2a\J  (crf  cos  (o;(t)))^  +  (e^  sin  {uj{t))f  (4-45) 

Oej2  =  2aed 

This  implies  that  even  under  J2  the  mean  relative  semi-major  axis  of  the  relative  for¬ 
mation  remains  equivalent  to  the  unperturbed  case.  Similarly,  the  radial  displacement 
is  nominally  found  as  Xd  =  5a,  and  assuming  mean  secular  rates,  this  value  remains 
unchanged  under  the  perturbation;  thus. 


Xdj2  =  5a  (4.46) 

The  in-track  displacement  is  effected  by  the  perturbation.  Substituting  state  expres¬ 
sions  into  the  yd  component  of  Eq.  4.43,  it  is  found  that 

ydj2  =  a59(t)  -|-  a5fl(t)  cosi  —  2aed  (sin  0  cos —  cos  6*  sin 0;^)  (4.47) 

Substituting  Eq.  4.44  into  Eq.  4.47,  using  the  dehntion  of  the  true  latitude  (6*  =  u+u) 
and  using  a  double-angle  trigonometric  identity  yields, 

ydj2  =  a-59o  +  a  {5uj  +  59)  /S.t  +  a[  (Jhlo  +  cos  i )  -I-  2aed  sin  {ujd.  —  9) 

V  ^  (4-48) 

=  a  (5^0+  +  a  ( (jci;  -|-  (5z>  -|-  cos  ij  At  +  2aed  sin  (ud  —  9) 

Immediately  it  is  observed  the  same  constant  offset  term  based  on  the  initial  argument 
of  latitude  and  ascending  node  is  present.  A  periodic  term  arises  that  induces  an 
oscillation  at  a  frequency  of  the  difference  between  the  deputy  perigee  and  chief  true 
latitude.  The  additional  term  to  the  in-track  offset  is  the  inclusion  of  the  perigee  and 
ascending  node  rates.  Secular  drift  in  track  is  manifest  in  this  condition,  which  implies 
that  59  =  5uj  +  59  =  —5VLcosi  will  limit  the  in-track  drift.  Einishing  the  in-plane 
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motion  with  the  perturbed  phasing  angle  and  substituting  the  perturbed  states 


tan  ^  I 

(  5qi  sin  9  —  6q2  cos  9 

\Sqi  cos  9  +  6q2  sin  9 

tan“^  I 

(  sin  {9  -  Ud)  \ 

\cos  {9  -Ud)) 

{9  —  Ud) 

(4.49) 


This  is  the  identical  expression  for  the  unperturbed  motion;  however,  the  argument 
of  perigee  term  now  has  a  secular  drift  from  the  J2  effect,  implying  that  the  frequency 
of  the  phasing  goes  as 

^J2  =  e  +  Ud  (4.50) 

Letting  9  for  the  circular  chief  be  expressed  as  n  +  ch,  the  in-plane  phasing  frequency 
is  found  as 

^J2  =  n-  6u  (4.51) 

where  dui  =  Ud  —  ojc- 

Focusing  in  the  cross-track  direction,  the  maximum  amplitude  is  found  as  just 
a  slight  modihcation  of  the  unperturbed  form  as 


Zniaxj2  =  a\  {6io)  +  sin  ioAt 


(4.52) 


Thus,  there  is  a  secular  drift  in  the  cross-track  direction.  This  can  be  eliminated  if 
the  ascending  node  rates  (<5^2)  are  matched.  If  not,  the  cross-track  magnitude  will 
increase  in  time  at  a  nearly  first  order  rate.  However,  for  a  chief  with  zero  inclination, 
the  cross-track  amplitude  will  remain  constant.  Finally,  the  cross-track  phasing  is 
found  as 


i!j2  =  tan’ 


I  SisinO  —  cos9smi6QAt\ 
cos  95i  +  sinO  siniSilAt  I 


(4.53) 
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Summarizing,  analytical  expressions  for  the  mean  response  of  the  ROEs  for  a 
circular  chief  are 


®eJ2  2oec( 

XdJ2  =  So- 

ydj2  =  o  (600  +  COS  i)  +  a  (Suj  +  6i>  +  6Q  cos  ij  At  +  2aed  sin  (ujd  —  0) 
/6j2  =  (^  ~  Old) 


(4.54) 


Zmaxj2  =  o^  (6io)  +  {  Sil  sin  ioAt 

'4)j2  =  tan 


I  I  SisinO  —  cos9  smiSQAt\ 
6i  cos  9  +  sin  9  sin  iStlAt  I 


These  equations  illustrate  the  importance  of  matching  the  drift  rates  of  the  chief  and 
deputy  ascending  nodes  and  arguments  of  perigee  for  formation  keeping. 


4-2. 1.3  Osculating  Relative  Orbit  Elements  from  the  J2  Effect  using  Or¬ 
bital  Element  Differences:  An  Analytical  Jacobian  Approach.  As  a  primarily  ana¬ 
lytical  exercise,  one  can  also  look  at  how  the  ROEs  will  change  as  a  function  of  time 
using  a  perturbation  approach  to  verify  the  previous  section’s  state  substitution.  Let 
the  solution  vector  X  equate  to  the  six  ROEs,  and  let  the  state  vector  5e  equal  to  the 
non-singular  orbital  element  differences  between  the  chief  and  deputy.  To  examine 
the  effect  of  J2  on  the  ROEs,  we  perform  the  following  Jacobian  operation 


(9X  _  (9X  d6e 
dt  d6e  dt 


(4.55) 


where 


dX 
06  e 


dae 

dae. 

dSa 

dSQ 

d-tp 

dip 

dSa 

■  ■  dsn 

(4.56) 


The  partials  in  Eq.  4.56  can  be  found  by  taking  partials  of  Eq.  4.43  with  respect 
to  the  orbital  element  differences.  For  this  exercise,  only  the  spatial  ROEs  will  be 
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examined,  while  the  phasing  angles  are  left  for  future  work.  The  mean  secular  J2 
rates  are  used.  The  time  rate  of  the  orbital  element  differences  is  expressed  as 


d6a 

dt 

0 

dse 

dt 

66 

dSi 

dt 

0 

dSqi 

dt 

qid 

dSq2 

dt 

q2d 

dsn 

_  dt  _ 

6Q 

The  resulting  matrix  math  needed  is  then 


dX 

d6e 


dae 

dSa 


d'ljj 

dSa 


0 

66 

dae 

dsn 

qid 

dtp 

q2d 

9(50  _ 

0 

6Q 

(4.57) 


(4.58) 


As  a  result,  only  the  second  through  fourth  and  sixth  columns  of  are  necessary. 
The  semi-major  axis  and  radial  displacement  are  found  as  functions  of  solely  a,  e^, 
and  6a]  thus,  partials  of  these  expressions  with  respect  to  any  other  variable  is  zero. 
Directly  Og  and  Xd  remain  constant-that  is. 


da^ 

~m 

dxd 

dt 


0 

0 


(4.59) 
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However,  examining  the  in-track  displacement  and  cross-track  amplification,  the  re¬ 
lationship  is  less  linear.  The  resulting  expressions  become 


dyd  9yd  ■ 

+ 


dyd 


dt 


Ozr 


dt 


859 

3z 

^^ma: 


-59  + 


d5qi 
dz. 


Qld  + 


dyd 

d5q2 


Q2d  + 


(9(5gi 


-Qld  + 


dZr, 


d5q2 


dyd 

d5Vt 

'-hd  + 


dZn 


(4.60) 


d5n 


'-5n 


Immediately  from  Eq.  4.43,  each  partial  of  yd  exists,  while  only  the  partial  of  Zmax 
with  respect  to  <514  is  non-zero.  The  gradient  of  yd  is  found  as 


dyd 

859 

dyd 

d5qi 

dyd 

d5q2 

dyd 

85n 


=  a 

=  —2a  sin  9 
=  2a  cos  9 
=  a  cos  i 


while  for  Zmax  we  find 


8z 


max 


d5VL 


a^^hlsin^  i 

^max 


Combining  the  above  results,  the  time  rate  for  yd  is  now 


(4.61) 


(4.62) 


dyd 

8t 


a59  —  2a  sin  9qid  +  2a  cos  9q2d  +  a  cos  i5fl 


=  a 


59  +  2  {q2d  cos  9  —  qid  sin  9)  +  512  cos  i 


(4.63) 


and  for  the  Zmax  time  rate 


dZr) 


a^5Vt  sin^  i 


8t 


5n 


And  if  we  express  512  =  5QAt,  Eq.  4.64 


d  Zrt 


8t 


sim  i{5^iyAt 


(4.64) 


(4.65) 
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Although  not  in  exact  agreement,  this  expression  shows  the  same  trend  in  the  Zmax 
parameter  when  compared  to  Eq.  4.54.  Equation  4.65  shows  a  linear  time-rate,  while 
Eq.  4.54  shows  a  relationship  on  the  order  of  The  two  expressions  appear  to 

have  time  rates  on  the  same  order  of  magnitude. 
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4-2.2  Relative  Orbit  Elements  for  the  Perturbed,  Circular  Chief  Developed  by 
the  Sehwieghart-Sedwiek  Differential  Equations  with  Arbitrary  Initial  Conditions. 
The  following  section  details  the  development  of  the  ROEs  from  solving  the  differential 
equations  and  setting  arbitrary  initial  conditions.  From  Eq.  3.12,  the  radial  and  in¬ 
track  motion  is  still  uncoupled  from  the  cross-track  as  in  the  HCW  model.  Placing 
the  in-plane  system  in  matrix  form  as  [36] 


X 

0 

0 

1 

0 

X 

y 

0 

0 

0 

1 

y 

X 

(5c^  —  2)n^ 

0 

0 

2nc 

X 

y 

0 

0 

—2nc 

0 

y 

(4.66) 


This  linear  system  can  be  solved  using  a  number  of  different  tools.  The  approach 
taken  by  this  author  is  the  matrix  exponential.  The  resulting  state  transition  matrix 
found  via  <h(t,to)  =  is 


$ 


SS 


B^—cf  cosh(ifi  At) 

0 

(  — S^+cr)  sinh(iCiA£) 

B—B  cosh(ifi  Ap 

K2 

(-K2P/2 

K2 

Bcr{K2At+^^^7^  sinh(ii:i  Ai)) 

1 

S(cosh(iCi  A£)  — 1) 

K2a  At+B"^  F/^^sinh(KiAt) 

K'2 

K2 

K'2 

—  K20-  smh.(KiAt) 

(-^2)3/2 

0 

cosh(iPiAt) 

-K2Bsmh{KiAt) 

(-i^2)3/2 

Ba{cos]i{KiAt)  —  l) 

0 

BiC2  sinh(Ki  Ap 

— cr+B^  cosh(iClA^) 

K2 

(-^2)3/2 

K2  j 

(4.67) 
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such  that  [x^y^x^yY  =  ^ss(t,to)[xo,yo,Xo,yoY ,  where  the  constants  are  defined  as 


At  =  {t-  to) 

Ki  =  Y-B‘^  +  a 
K2  =  B‘^  -a  =  -Kl 
B  =  2nc 


a  =  (5c^  —  2)n^ 


3J2RI 

8r2 


(1  +  3  cos(2i)) 


(4.68) 


There  now  exists  a  closed  form  solution  of  the  SS  differential  equations  regardless  of 
initial  conditions. 


4-2.2. 1  Deriving  the  In-Plane  Relative  Orbit  Elements.  Focusing  on 
the  in-plane  motion,  the  state  can  be  expressed  linearly  at  any  time  as  X{t)  = 
^{t,to)Xo.  From  Eq.  4.16,  the  initial  Cartesian  conditions  were  found  as  functions  of 
the  initial  ROEs  by  the  following  relation 

1  0 

sin  /3o  0  1 

I  sin  /3o  0  0 

n  cos  Po  0 


Allowing  the  product  of  time-varying  matrix  in  Eq.  4.69  with  the  STM  in  Eq.  4.67  to 
be  expressed  as  J(t),  the  instantaneous  ROEs  are  found  as  R(t)  =  A  (J(t)i?jo),  where 

r  1  ^ 

A  is  the  transformation  from  Cartesian  to  ROE  space,  and  Rio  =  Oeo  Xdo  ydo 
the  initial  in-plane  spatial  ROEs.  Now  knowing  the  necessary  operations,  the  only 
remaining  issue  is  the  algebraic  derivation.  For  simplicity,  it  is  desired  to  express  the 
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in-plane  relative  position  in  the  following  form 


X  =  Cl  +  Cxi  cosh  6  +  Cx2  sinh  6 
y  =  C2  +  Cy\  cosh  9  +  (y2  sinh  9  +  (pAt 


(4.70) 


The  expression  in  Eq.  4.70  is  the  X  =  J{t)RiQ  function.  Performing  the  matrix 
multiplication  and  grouping  like  terms,  the  constants  in  Eq.  4.70  are  found  as 


Cl  =Xdo{^ 

R-2 


2K2 

C2  =  Vdo  +  fleo  (  sin/3o 
/  3Bn 


3Bn\  /  5  cos /do  cos /do 


Cxi  Xdo 
Cx2  = 


\2K2 

tteon  sin  /do 


a 

1^2 


Bn  sin  /do 
2X0 


BcosPo  acos/3o] 

+  fleO  I - 77 - K 


K2 


2Ko 


Cyl  - 


2^^2 

aeoBnsin  /do 


(4.71) 


2X0 


Cy2  —  X  do 

=  XdO 
9  =  KiAt 


3B^n 


Ba 


2{-K2fB 

Ba  3na 


-K2YC 

I"  fleO 


+  fleO 


B‘^  cos  /do  Ba  cos  /do 


+ 


{-K2YC  2(-X2)3/2 


a  COS /do  ^  Ba  cos  Po 


K2 


2K2 


Also  from  Eq.  4.70,  time  derivatives  provide  the  relative  velocities  to  be 


X  =  9  (Cxi  sinh  9  +  Cx2  cosh  9) 

y  =  9  (Cj/i  sinh  9  +  (y2  cosh  9)  +  Lp  (4.72) 


where  9  =  Ki. 

Having  condensed  forms  of  the  relative  trajectory  expressions,  the  in-plane  por¬ 
tion  of  the  A  transformation  is  now  possible.  Beginning  with  the  semi-major  axis  Oe, 
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define  the  following  two  variables 


X 

del 

n 


=  -  {Cx2  cosh  e  +  Cxi  Sinh  9) 
n 

ae2  =  3x  +  2  — 
n 


=  3  (Cl  +  Cxi  cosh  9  +  Cx2  sinh  9)  -\ - (Cj/2  cosh  9  +  C?;i  sinh  9) 


such  that  the  final  expression  for  becomes  by  simple  substitution 


(4.73) 


de 


2\/ d^i  +  d 


^e2 


(4.74) 


which  is  seen  as  the  combination  of  hyperbolic  trignometric  functions  at  a  frequency 

oil 

The  in-plane  displacements  along  the  radial  and  in-track  directions  can  be  writ¬ 
ten  as  a  linear  combination.  Using  the  expressions  from  Eq.  2.59 


Xd  =  4:X  - 

n 

2x 

yd  =  y - 

n 


(4.75) 


Substitution  of  Eq.  4.70  and  Eq.  4.72  and  grouping  terms  according  to  the  trigono¬ 
metric  terms,  it  is  observed  that 


Xd 

yd 


4Cxl  +  ^  4Cx2  + 

2ec,:,2 

n 


Cwl 


Cy2  - 


2gC»i 

n 

2gC^i 

n 


cosh0 

sinh6' 


4Ci 
C2  + 


(4.76) 


The  in-plane  phasing  is  hnally  found  as 


f3  =  tan 


-1 


X 


3nx  +  2y 


(4.77) 
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and  using  the  previously  defined  variables  Ogi  and  ae2  from  Eq.  4.73,  we  find 

=  tan“^  (4.78) 

V«e2/ 

4- 2. 2. 2  Out-of-Plane  Relative  Orbit  Elements.  Schweighart  and  Sed- 
wick  [37]  further  defined  the  cross-track  motion  for  the  circular  satellite  under  the  J2 
perturbation.  The  cross-track  motion  is  described  in  [37]  as 

z{t)  =  A{t)  sin  {B{t)t  —  C{t))  (4.79) 


where  Schweighart  denotes  A{t)  as  the  time- varying  cross-track  magnitude,  B{t)  as 
the  orbital  frequency  (which  is  modified  under  J2),  and  C(t)  as  a  phasing  angle. 
Expressed  as  functions,  these  parameters  are 


A{t)  =  rrj 
B{t)  =  nk 

(4.80) 

k  =  c  +  I  cos^  i 

^  3J2RI 

2r2 

An  immediate  comparison  can  be  made  to  the  ROE  expression  for  the  cross-track 
motion 

z(t)  =  ZmaxAnfj  (4.81) 

such  that  it  is  evident  that 


-ipit)  =  B(t)t  —  C(t) 


(4.82) 


The  time-varying  magnitude  A{t)  is  found  as  the  product  of  the  circular  reference 
orbit’s  radius  and  the  spherical  angle  [r])  resulting  from  the  combination  of  differences 
between  the  ascending  nodes  and  inclinations  of  the  chief  and  deputy.  For  small 
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angles,  this  has  already  been  found  earlier  in  Eq.  4.43  as  Zmaxh'  and  is  equal  to 


rj  =  y  (Ai)^  +  (Aflsin^i)" 


(4.83) 


and  develops  according  to 


7]  =  rjo  cos  7o  sec  7 


(4.84) 


where  7  is  the  cross-track  phasing  angle,  expressed  for  small  angles  as  [37] 


7  =  cof 


-1  ( 

\  Af2  sin  i 


(4.85) 


Using  the  proposed  parameterization  in  Eq.  4.82,  and  assuming  to  =  0,  the  angle  7 
is  equal  to  and  opposite  tjj  at  epoch.  That  is 


^0  =  -70 


(4.86) 


which  is  given  by  Schweighart  and  Sedwick  as 


.  -I  f  zo  1 


(4.87) 


where  b  =  I  sin^  i.  This  immediately  implies  from  Eq.  4. 


^0  =  tan  ^  (—  ^ 

\2:o  n[K  —  b) 


(4.88) 


The  cross-track  phasing  is  then  given  as 


tjj  =  nkt  +  Tpo 


(4.89) 


Turning  the  attention  to  the  cross-track  amplihcation,  the  expression  for  A{t)  is  given 


A{t)  =  rr] 


(4.90) 


=  rrjQ  cos  7o  sec  7 
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Utilizing  Zmax  =  ^(t)  and  making  the  substitutions  Zmaxo  =  and  'ipo  =  —70,  we 
arrive  at 

A(t)  =  Zmaxo  cos  'ipo  sec  7  (4.91) 

where  7  is  dehned  in  [37]  as 

7  =  tan“^  {nbt  +  tan  70)  (4.92) 

Finally,  from  Eq.  4.82,  the  cross-track  amplihcation  is  found  as 

Zmax{t)  =  ZmaxO  COS  SCC  7  (4.93) 

such  that  the  cross-track  motion  is  then  dehned  as 

z{t)  =  Zmaxo  COS  V’o  SCC  7  siu  ^7(t)  (4.94) 

Utilizing  the  trigonometric  identity 

sec  (arctan(a;))  = - 7 

cos  (arctan(a;)) 

=  ^  (4.95) 

=  VT+Ir^ 

and  substituting  in  Eq.  4.94,  the  hnal  form  of  the  cross-track  motion  is  then  found 
as 

z{t)  =  Zmaxit)  (4.96) 

where 

Zmax{t)  =  ZraaxO  COSIpox/ 1  +  {but  -  tanipof  (4.97) 
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4 -2. 2. 3  Summary  of  the  Relative  Orbit  Elements  for  the  J2  Perturbed, 
Circular  Chief  for  Arbitrary  Initial  Conditions.  The  following  provides  a  summary 
of  the  ROEs  for  the  circular  chief  perturbed  by  J2- 


a-e  —  ‘2y  ah  +  ah 

Xd  _  4:Cxi  +  4Cx2  +  cosh  9  ^  4Ci 

Vd]  42-%J  [sinh^J  [c2  +  cpAt 

(3  =  tan“^ 

\Oe2/ 

^max  {t)  =  Zmaxo  cosfjQsJ  1  +  {but  -  tanipof 
fj  =nkt  +  fjo 


(4.98) 


These  expressions  have  been  developed  for  arbitrary  initial  conditions  and  are  not 
dependent  on  velocity  states  initialized  to  bound  the  relative  orbit. 


4-2.3  Numerical  Examples  and  Deviations  from  the  Clohessy-Wiltshire  As¬ 
sumptions.  As  the  primary  assumption  in  the  HCW  realization  is  a  circular  chief, 
examining  the  effects  on  the  HCW  ROEs  from  J2  is  very  appropriate.  The  stationary 
relative  orbits  are  found  using  the  conditions  detailed  in  Appendix  F.  The  following 
numerical  examples  will  demonstrate  the  deviation  of  the  model  from  HCW.  Initial 
conditions  are  given  in  the  form  Rq  =  Oeo  Xdo  Vdo  Po  z^axo  f’o  units  of 

(km, rad)  where  appropriate,  and  the  chief  orbital  elements  are  given  in  the  form 
Gco  =  a  e  i  u  Vt  Mq  with  units  of  (km,rad)  as  well. 

4.2.3. 1  Relative  Orbit  Elements  at  for  a  Low  Altitude  Chief  with  Low 
Inclination.  Given  the  following  initial  conditions 

Ro  =  [ll.89  0.01  0.05  0  1  0.5 

:  ^  (4.99) 

eco  =  7000  0  0.1745  0  0  0 
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Figure  4.10  provides  a  visualization  of  the  in-plane  ROEs.  Initial  conditions 
for  bounded  motion  have  been  used.  The  expressions  for  these  initial  conditions  are 
derived  in  App.  F.  Of  immediate  note  is  these  initial  conditions  supplied  to  the  HCW 
model  yield  a  drifting  trajectory  (evident  in  There  are  variations  on  the  order  of 
10^  meters  in  the  expressions  for  OejXrf,  and  Ud-,  however,  the  in-plane  phasing  differs 
only  slightly  from  the  expected  motion. 


Figure  4.10:  In-Plane  Osculating  Relative  Orbit  Elements  for  J2  Perturbed,  Low  Al¬ 
titude  Circular  Chief 


Figure  4.11  provides  a  visualization  of  the  out-of-plane  ROEs.  The  cross-track 
amplihcation  decreases  nearly  linear,  but  at  low,  near  negligible  rate.  This  again  is 
due  to  the  bounding  condition.  Differential  drift  rates  will  induce  perturbations  to 
the  orbital  angular  momenta  of  both  the  chief  and  deputy  and  modify  the  action 
of  the  cross-track  motion.  The  short  term  effect  of  the  cross-track  variation  is  near 
negligible,  but  long  term  projections  will  be  error  prone. 


Figure  4.11:  Out-of-Plane  Osculating  Relative  Orbit  Elements  for  J2  Perturbed,  Low 
Altitude  Circular  Chief 
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4-3  Applications  of  Osculating  Relative  Orbit  Elements 


4-3.1  Guidance  and  Navigation. 


4. 3. 1.1  Derivation.  Having  analytical  expressions  for  the  osculating 
ROEs,  one  can  examine  the  effect  of  an  impulsive  burn  in  any  of  the  directions  in  the 
orbital  plane  of  the  chief.  The  derivation  begins  by  assuming  that  a  single,  impulsive 
burn  is  possible  and  can  be  oriented  in  such  a  way  that  the  burn  acts  entirely  in  the 
radial,  in-track,  and/or  cross-track  directions.  The  derivation  follows  similarly  to  that 
in  [33],  and  focuses  primarily  on  the  idea  of  a  single  burn.  Assuming  knowledge  of 
the  ROE  state  at  the  time  the  burn  {th)  occurs,  and  letting  the  subscript  pb  indicate 
the  desired  post-burn  value,  the  expressions  can  be  written  as 


®epi, 

^etb 

AOe 

^dpi, 

^dbb 

Axd 

Vdpi 

Vdtb 

+ 

Apd 

fdph 

fdtb 

Aft 

^maXpi, 

\z 

^^max 

'^ph 

i^tb 

Alp 

(4.100) 


where  the  A  values  indicate  the  contribution  to  the  ROE  provided  by  the  impulse.  If 
the  relative  velocity  is  known  at  t  =  tb,  then  an  instantaneous  burn  would  yield  the 
following  expressions  for  the  A  components 


1 

<1 

<1 

^Db  +  ^ 

Apd 

ytb-^ 

A(4 

'  (3..»,AAd 

\z 

^^max 

\/(^) 

Alp 

tan-> 

(4.101) 


where  the  A  values  are  expressed  as 


Ax 

X  +  AI4 

Ay 

= 

y  + 

Az 

i  + AI4 

(4.102) 


One  can  then  make  the  substitution  using  the  Cartesian  to  ROE  conversion  such  that 


ae_  cos  /3_ 

- 

ytb  =  Oe-  sin  /3_  +  ya-  ^'^3) 

Ztb  =  Zmax-  sin  tjj. 

where  the  minus  subscript  indicates  the  state  immediately  before  the  burn.  Substitu¬ 
tion  into  Eq.  4.101  yields  dependence  entirely  on  the  ROE  states.  If  the  expressions 
for  the  A  values  in  Eq.  4.101  represent  desired  maneuvers,  Eq.  4.101  represents  a 
system  of  highly  non-linear  equations  that  can  be  numerically  solved  to  determine 
approximately  the  necessary  velocity  change  needed  to  alter  the  relative  trajectory. 


4 .3. 1.2  Mathematical  Concerns.  It  is  quite  signihcant  to  note  that 
the  nonlinear  equation  set  in  Eq.  4.101  suffers  from  constraint  issues.  For  example, 
if  a  desired  Og  is  selected,  there  are  two  independent  variables  {Ax,  Ay)  to  perform 
the  maneuver.  The  out  of  plane  motion  is  such  that  only  a  Az  will  alter  the  cross¬ 
track  trajectory.  If  Aog  is  desired,  an  inhnite  amount  of  solutions  could  exist  due 
to  the  two  free  variables  in  the  expression.  This  implies  that  if  a  Aog  is  selected, 
then  one  of  either  Axd,  Aya,  Ap  should  be  selected  additionally.  Selecting  either  Axd 
or  Ayd  will  yield  a  possible  solution,  as  does  the  selection  of  both  parameters.  It  is 
nearly  mathematically  impossible  to  select  both  Ap  and  Azmax  unless  selected  at  a 
location  in  which  ztb  =  0.  Based  on  these  constraints,  an  arbitrary  maneuver  is  not 
mathematically  possible.  There  are  only  3  independent  variables  with  6  dependent 
variables.  One  can  only  reach  a  subspace  of  the  desired  maneuver  for  any  given 


relative  state  at  time  tj,.  The  coupled  nature  of  the  variables  is  investigated  in  the 
next  section. 


4 .3. 1.3  Numerical  Examples.  It  is  instructive  here  to  examine  the  ef¬ 
fects  of  impulsive  burns  on  the  analytical  expressions.  The  entire  COE  set  of  the  chief 
is  not  needed  as  the  only  value  needed  from  the  derivation  is  the  mean  motion.  The 
ROEs  have  been  treated  as  parameters  that  can  be  acquired  via  relative  navigation 
data,  and  the  only  needed  values  are  those  prior  to  the  burn. 

For  a  chief  of  any  eccentricity  representing  a  closed  orbit  with  a  semi-major  axis 
of  7000  km,  the  following  burn  time  ROEs  are  given  as 


(Xg— 

0.01  km 

Xd- 

0.01  km 

Vd- 

0.02  km 

(3- 

TT  rad 

^max— 

0.01  km 

_ 

0.1  rad 

(4.104) 


Allowing  the  impulses  in  the  radial  and  in-track  directions  to  vary  between  - 
1  and  1  m/s.  Fig.  4.12  provides  a  visualization  of  the  coupling  effect  on  the  ROE 
response. 

Immediately  evident  is  the  linear  response  in  the  values  for  Xd  and  i/d.  Also,  a 
linear  response  is  noted  in  the  response  for  z^ax-  The  most  non-linear  is  the  response 
of  the  semi-major  axis,  which  forms  an  open  shape  resembling  a  paraboloid  that  grows 
withont  bonnd.  Resulting  spatial  responses  regardless  of  initial  conditions  take  the 
same  form  as  a  conseqnence  of  the  deterministic  derivation. 
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Figure  4.12:  Coupling  Among  the  Spatial  ROEs  with  Respect  to  Three-Dimensional 
Impulsive  Velocity  Burns 
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V.  Conclusions  and  Future  Work 


5.1  A  Dialogue  on  the  Realism  of  Osculating  Relative  Orbit  Elements 

Within  this  thesis,  a  method  has  been  shown  to  parameterize  the  relative  tra¬ 
jectory  of  two  satellites  in  close  proximity  using  geometric  insight.  The  effects  of  the 
chief  eccentricity  and  the  J2  perturbation  on  the  circular  chief  have  been  investigated. 
For  the  latter,  mean  J2  effects  have  been  found  to  serve  as  a  main  variable  of  drift 
dependence.  Perturbations  in  the  relative  orbits  will  only  further  compound  in  the 
relative  dynamics.  If  attempting  to  attain  geometrical  insight,  full  inclusion  of  the 
chief  eccentricity  forces  the  model  into  the  domain  of  mathematical  abstraction;  how¬ 
ever,  the  relative  orbit  element  (ROE)  parameterization  still  serves  as  a  methodology 
to  determine  the  instantaneous  magnitude  and  phasing  of  the  oscillatory  motion,  in 
addition  to  a  time- variant  offset. 

Relative  motion  is  an  abstract  concept,  and  quite  difficult  to  visualize  without 
in-depth  study.  Oscillating  relative  trajectories  from  differential  orbits  and  energy 
exchange  are  difficult  concepts  to  grasp.  More  so,  even  the  time-varying,  rotating 
local-vertical,  local- horizontal  (LVLH)  orbital  frame  in  which  the  relative  dynamics 
occur  is  an  abstract  notion.  These  factors  compound  and  convolute  mission  critical 
objectives  including  relative  navigation,  guidance  and  control,  and  docking  maneu¬ 
vers.  The  idea  is  to  obtain  a  geometrical  parameterization  that  allows  an  operator 
to  quickly  visualize  the  relative  motion  and  further  the  mission  of  close  proximity 
operations.  However,  the  osculating  eccentric  ROEs  fail  to  meet  this  objective  in  this 
study;  the  mathematical  notion  of  a  locus  of  points  satisfying  a  2  x  1  ellipse  does  not 
match  criteria  of  operational  efficacy.  For  LEO  missions,  the  use  of  a  J2  perturbed  set 
will  lead  to  higher  accuracy  in  navigation;  in  contrast,  the  physical  near  impossibility 
of  a  true  circular  orbit  may  lead  to  eccentricity  effects  dominating  the  perturbations. 

In  this  study,  the  mathematics  have  been  done  to  fully  show  how  the  ROE 
parameterization  grows  over  time  with  state  substitution  relaxing  the  unperturbed 
and  the  circular  chief  assumption.  The  motion  was  found  to  be  fully  mapped  by 
these  analytical  expressions.  Moreover,  closed-form  expressions  for  the  relative  orbit 
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elements  to  include  these  effects  now  exist.  Further  work  is  still  needed  to  develop  a 
fully  intuitive  and  geometrically  insightful  parameterization  to  visualize  and  control 
the  relative  motion  of  bodies  in  close  proximity. 

5.2  Recommendations  for  Future  Work 

Throughout  the  course  of  this  study,  several  alternate  routes  have  been  found 
that  may  be  followed  in  the  future.  These  may  help  further  the  use  of  the  ROEs  to 
hnally  develop  a  parameterization  encompassing  eccentricity  and  the  J2  perturbation. 

5.2.1  Velocity  Independence.  The  expressions  for  the  ROEs  include  substi¬ 
tution  of  the  relative  velocity  states.  If  one  is  comparing  ROE  states  across  various 
perturbed  environments,  the  rotation  of  the  LVLH  frame  induces  different  angular 
velocities.  These  angular  velocities  modify  the  velocity  states  signihcantly.  The  non¬ 
linear  operations  on  the  relative  velocity  terms  compound  differences.  A  proposed 
route  is  now  to  investigate  the  analytical  HOW  expressions  to  hnd  ROE  expressions 
that  are  simply  functions  of  the  relative  position  terms.  This  has  been  done  using 
orbital  element  differences.  The  estimation  of  inertial  positions  and  velocities  with 
conversions  to  two  sets  of  orbital  elements  is  possible;  however,  linearizations  in  the 
model  limit  the  range  of  applicability  of  this  method.  Therefore,  it  is  highly  desired  to 
identify  relative  orbit  elements  as  functions  of  position.  Furthermore,  it  is  recommend 
that  a  parameterization  be  made  that  is  based  on  a  model  developed  in  inertial  space. 

5.2.2  Mathematical  Inspection  of  the  Yamanaka- Ankers en  Topology.  Simi¬ 
lar  to  re-parameterizing  the  ROEs  based  on  relative  position  terms,  the  Yamanaka- 
Ankersen  model  needs  to  be  investigated  and  placed  into  a  form  to  understand  the 
resulting  topology.  The  uncoupled  out-of-plane  and  in-plane  motion  should  be  ana¬ 
lyzed  separately  and  later  superimposed.  This  was  somewhat  attempted  during  this 
study  with  little  gain.  Placing  the  closed  form  solution  into  an  expression  that  is  a 
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function  of  the  eccentricity  of  the  chief  in  addition  to  initial  conditions  could  motivate 
topological  intuition. 

5.2.3  Higher- Order  Terms  in  the  Virtual  Chief  Model.  Propagation  of  the 
relative  motion  using  a  virtual  chief  approach  could  yield  intuitive  results  if  the  accu¬ 
racy  of  the  model  were  increased.  A  linear  state  transition  matrix  was  developed  by 
Johnson  [27]  but  was  found  to  be  error  prone.  The  idea  of  retaining  higher  order  terms 
in  the  close  proximity  linearization  of  the  HCW  model  would  undoubtedly  increase 
the  accuracy  of  the  dual  propagation. 

5.2.4  Drag  Effects  on  the  Relative  Orbit  Elements.  Inherent  in  the  assump¬ 
tion  in  the  ROE  development  is  the  idea  of  near  identical  ballistic  coefficients  yielding 
identical  drag  effects.  An  interesting  route  of  study  is  to  examine  analytically  the 
effects  on  the  ROEs  including  the  effects  of  atmospheric  drag.  As  certain  ROEs  have 
been  found  to  be  functions  of  orbital  element  differences,  the  induced  energy  decay 
of  the  orbits  would  also  effect  the  ROE  parameters. 

5.2.5  Perturbation  Methods.  Rather  than  performing  full  state  substitution 
of  closed-form  solutions  into  the  ROE  parameterization,  a  worthwhile  study  would 
be  to  examine  the  effect  of  the  ROEs  using  the  following  method.  Let  R  denote  the 
ROE  states,  and  X  denote  the  relative  Cartesian  state.  The  differential  states  can 
then  be  examined  by 

JR=— JX  (5.1) 
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Appendix  A.  State  Transition  Matrix  Properties 

The  state  transition  matrix  (STM)  is  a  deterministic  method  of  propagating  a  state 
at  time  to  to  any  time  t.  The  following  identities  and  concepts  are  presented  as 
fnndamental  and  withont  proof. 

Typically,  the  propagation  resembles 

X(t)  =  $(t,to)X(to)  (A.l) 

where  X  is  the  state  vector  and  <h  is  the  STM.  For  a  linear  time-invariant  system  of 
the  form 

X  =  Ax  +  Bu  (A. 2) 

where  A  is  the  plane  matrix,  B  is  the  control  matrix,  and  u  is  the  inpnt  vector,  the 
solntion  is  given,  withont  replication  of  proof,  as 

X{t)  =  +  f  e^^^-^Mu{T)dT  (A.3) 

Jto 

with  r  as  a  dnmmy  variable  of  integration.  The  integral  in  Eq.  A.3  is  often  termed 
the  convolntion  integral.  The  matrix  exponential  in  Eq.  A.3  is  often  represented  as 

e^('-*o)  =  (A.4) 


For  any  order  of  linearity,  the  STM  must  satisfy  the  following  properties 


1.  <h(t2, to)  —  *h(t2, ti)<h(ti, to) 

2.  <h(t2, ti)  =  <h“^(ti, t2) 

3.  4>(ti,ti)  =  1 

Also  as  a  direct  consequence  of  Eq.  A.4,  the  STM  must  also  satisfy 


<9$(t,to) 

dt 


A<h(t,to) 


(A.5) 
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Appendix  B.  Derivation  of  the  Eeeentrie  Oseulating  Relative  Orbit 

Elements 


The  following  appendix  provides  a  more  detailed  derivation  of  the  eccentric  ROEs 

iT 


using  state  substitution.  Dehning  the  in-plane  state  vector 
known  from  Sec.  4.1  that 

X  =  B(t)Ajpo 

iT 


X  y  X  y 


as  X,  it  is 


(B.l) 


where  Ajpo  =  |a,(, 
dehned  in  Sec.  4.1  as 


,  and  the  4x3  matrix  B(t)  is  described  in  variables 
-Bn  •  •  •  -Bis 


B(()  = 


(B,2) 


B4I  .  .  .  B43 

Combining  the  above  expressions,  expressions  for  the  states  are  found  carrying  out 
the  multiplication  as 

X  =  BiiOeO  +  Bi2XdQ  +  Bi^ydO 

y  =  B2iae0  +  B22XdO  +  B23yd0 

X  =  BsiOeo  +  B^2Xdo  +  B^^ydo 
y  =  B^itteo  +  B^2Xdo  +  B^^ydo 


(B.3) 


Now,  making  use  of  our  state  substitution  method,  the  ROE  expressions  can  now  be 
found. 


B.l  Relative  Semi-major  Axis,  Oe 

From  the  previously  described  relation,  Oe  is  found  as 


Ojq 


(B.4) 
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Substituting  in  the  expressions  for  x,  x,  and?/,  we  separate  the  terms  in  parenthesis  as 


BsitteO  +  B^2XdO  +  B^sUdO 


X 

n )  n 

o  I  ‘^y\  —  f  o  u  I  “^Bai  ^ 

3a;  H - 1  —  I  oBn  H - I  Oeo 

n  /  \  XI  J 


3i?i2 


2^42  \ 
n  ) 


■?R  ^2i?43\ 

3^13  H - <^e0 

n  / 


Squaring  the  two  expressions  in  Eq.  B.5 

VdO^h  I 


^2  r2 
®e0^31 


^2  r>2 

Xd0X^32 


2ae03;  £(0-831  ^32  ,  “^deOUdoBsiB^^  ‘^XdoydoB^2B3,‘i 


3a; = 


-'■eO 


^11 


4S|i 


Xdo  (  9-8^2 


4B 


2  /  f)p2  I  4-843 

J/do  1  9-8i3  + 


(^eO^dO  ISB11B12  ■ 


12B11B41 
v? 

12-Bi3i?43\ 

) 

8-B41-B42  12i?i2i?41  12_Bii_B42 


.I2  12^12-842  \ 

V?  V?  J 


rB 

+ 

n 

+ 

n 

8B4ii?43 

+ 

12i?i3i?4i 

+ 

12Biii?43 

Tn? 

n 

n 

8B42B43 

+ 

12^13^42 

+ 

12^12-643 

rB 

n 

n 

(B.5) 


(B.6) 


Now,  collecting  like  terms  in  Eq.  B.6  with  respect  to  the  initial  ROEs,  we  can  End 
the  a  coefficients  described  in  Sec.  4.1  and  detailed  in  Appendix  C.  The  resulting 
groups  need  to  be  multiplied  by  a  factor  of  4  to  compensate  for  the  multiplication  by 
2  in  Eq.  B.4.  For  example. 


0^1 


V  n2 


+  9Ri\+4^  +  12  \  ) 

J 

2  R?i  RiiR4i\ 

+  36^4^  +  16^  +  48 

J 


and  the  process  would  follow  by  inspection. 


(B.7) 


B.2  Radial  Displacement,  Xd 

The  radial  displacement  is  expressed  as 


Xd  =  + 


2y 

n 


(B.8) 
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This  can  be  found  through  simple  substitution  as 


/I  \  _  A  (U  ID  ID  'ii^  (-^41060  +  Bi2XdO  +  -B43l/(io) 

4x  H - —  4  \Biiaeo  +  Bi2Xdo  +  -Bisl/do)  H - 

n  n 

=  ( 4Bn  +  +  ( 4B.2  +  +  f 4B,3  +  !,«  (B.9) 

\  n  J  \  n  /  \  n  J 

=  (Jiaeo  +  o'2a;rfo  + 

5.5  In- Track  Displacement,  yd 

Very  similar  to  the  radial  displacement,  the  in-track  displacement  is  expressed 
as 

2t 

yd  =  y--  (B.io) 

n 

and  simple  substitution  yields 
2x 

yd  =  y - 

n 

/  D  ,  D  I  D  ^  ^  (-Bsi^eO  +  B^2XdO  +  B^^ydo) 

—  {B2iaeQ  +  B22XdO  +  -D23l/do) - 

n 

—  (  O  ^  2i?32^  ^  ^  2i?33^  IB  111 

—  I  -021 - a-eo  +  I  B22 - Xdo  +  I  B2S - ydo  ^  ' 

\  n  J  \  n  J  \  n  J 

=  SiOeO  +  '^2XdO  +  ^sydO 


Interestingly,  the  two  ROEs  that  are  linear  combinations  remain  in  a  similar  form  as 


Xd 

0-1 

0-2 

0-3 

yd 

El 

E2 

E3_ 

fleO 

^dO 

VdO 


(B.12) 


B.4  In- Plane  Phasing  ,  (3 

The  in-plane  phasing  is  also  found  by  simple  state  substitution.  Expressing  (3 
as 

f3  =  tan“^ 
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3>nx  +  y 


(B,13) 


Having  already  found  expressions  similar  to  the  numerator  and  denominator  while 
hnding  the  Oe  expression,  the  first  of  Eq.  B.5  is  multiplied  by  n,  and  the  second  of 
Eq.  B.5  is  divided  by  n  such  that  the  resulting  expression  is 


f3  =  tan 


-1 


BsitteO  +  By,2XdO  +  -Bssl/dO 


(3ni?ii  +  2H41)  Oeo  +  (3ni?i2  +  2H42)  XdQ  +  (3ni?i3  +  2H43)  ydo 


(B.14) 


B.5  Out-of-Plane  ROEs 

The  out-of-plane  ROEs  were  derived  in  Sec.  4.1  and  followed  a  rather  simple 
development. 
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Appendix  C.  Time-  Varying  Parameters  of  the  Yamanaka-Ankersen 

Derived  Relative  Orbit  Elements 

The  following  describes  analytically  the  time- varying  coefficients  of  the  relative  orbit 
elements  for  an  eccentric  chief.  The  coefficients  are  presented  assuming  an  initial  true 
anomaly  at  perigee  and  epoch  time  of  zero. 


C.  1  Matrices 


C.1.1  Members  of  the  A  Matrix. 


An  = 


( — 1  —  e)  ( — 3c  (1  -f-  e)  -T  ( — 1  -1-  -|-  3  (1  -t-  e))  (2  —  3eJ s)  A 


1  -e2 


Aio  —  0 


^13  — 


Au  — 


(-c  +  2e) 

(1  —  e^)  D 

C  (c  (e  -f  2  -|-  2e)  C  —  (1  +  e)^  (2  —  3e  J s)  C) 


(1  —  eA 


A21  —  0 


(C.l) 


A22  —  (1  +  e)  C 

C  (2  -  e  -  e^)  C  -  c  (c  -  2e)  AC 


A23  — 


^24  — 


(1  —  eA  k'^ 

C(-^^^  +  (e  +  2  +  2e)  ACs) 
(1  —  eA  D 
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^31  —  —  - 2  ( — ^  ( — 3(1  + c)  {c"y  +  c^')  +  ( —  l  +  e^  +  3(l  +  e))  ("7(2  —  3eJ s)  —  3e^'(t/s  +  sC^) ) ) 


A32  =  0 
^33  =  — 
A34  =  — 

A41  = 


(— c  +  2e)  (js  +  'I' s')  ( 


(1  —  e^) 


1 


(1  —  e^)  k^ 
1 


e  +  2  +  2e)  (cj  +  c''!')  —  (1  +  e)^  (7  (2  —  3eJs)  —  3e'I'  ( Js'  +  sC^)) 


— 1  —  e)  ^ — 3  (1  +  e)  ((2c  —  e)  ^'  +  7^'S)  +  ( — 1  +  +  3  (1  +  e))  ^3^'  (1  —  2eJ s)  H — ^ 


1  — 

^42  =  (1  +  e)7 

((2  —  e  —  e^)  7  +  (c  —  2e)  (— C7A  +  2'1's)) 


^43  =  — 


2144  = 


c 


(1  —  e^)  k"^ 

(e  +  2  +  2e)  ((2c  -e)'i>  +  -fXs)  -  (1  +  ef  (^3^'  (1  -  2eJs)  + 


(1  —  e^)  k^ 

where  the  variables  A,  7,  (,  and  ^1/  are 


V 


3jJ\ 

e ) 


(C.2) 


7  =  k‘^esC, 

d/  =  k^/c 
C  =  i/p 
A  =  l  +  C 


(C.3) 
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C.1.2  Members  of  the  B  Matrix. 


Bu  =  A.,  sin  A  -  co=  A  + 

-Bis  =  ^12 

D  .  .  ^  Aaicos/^o  A23nsin/3o 

B21  =  A22sm(3o - - - h  A24ncos(3o  H - - - 


3nA 


24 


B22  —  A21 - - - 

-B23  =  ^22 

D  _  4  ■  n  ^31COS/3o  ,  .  .  ,  Aggnsin/^o 

-B31  —  ^32  sin  Po - - - h  ^34^  cos  Po  H - - - 


-B32  —  ^31  ~ 


3nA, 


34 


-B33  —  ^32 

D  A  ■  o  ^4icos/3o  A43nsin/3o 

B41  =  /I42  sin  /:/o - - h  ^44^  cos  Po  H - - 


-B42  —  ^41  ~ 


3nA 


14 


-B43  —  ^42 


(C.4) 


C.1.3  Members  of  the  C  Matrix. 


Cu  = 


C12  — 


PePAe 
SAe 

k^PoPePAe 


(cAePo  —  esoSAe) 


C21  —  Pok 


2  ( CAoSe^  —  Pe-SAe^ 


V 


PePAe 


J 


es  k"^  ( ^^^ke  +  ^^AeSe 
\  PePAe 


C22  —  — 


f  CAepj  +  egAgsA 


Po  V 


PePA6» 


y 


(C.5) 


C.1.4  Members  of  the  D  Matrix. 


Dll  =  Cn  sin  i/jq  +  Ci2n  cos  i/jq 
D21  =  C21  sin  fjQ  +  C22n  cos  i/jq 


(C.6) 
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C.2  Time  Varying  Coefficients 


C.2.1  a  Coefficients. 

ABl  16^2  485ii54 


’ll 


’12 


’l3 


+ 

+ 

n 

+ 

I6-8I2 

+ 

48-812-842 

n 

45|3 

+ 

IQBl, 

+ 

48-813-843 

n 

D  D  8-B32-B33  32542-B43  48-Bl,3i?42  48-812-84 

0:4  =  72-8i2-8i3  + 


as  =  725n5i2  + 
ae  =  72-8ii-8i3  + 


+ 

+ 

n 

+ 

n 

8-831-832 

+ 

32-841-842 

+ 

48-812-841 

+ 

48-811-842 

n 

n 

8-831-833 

+ 

32-841-843 

+ 

48-813-841 

+ 

48-811-843 

(C.7) 


n 


n 


(T  Coefficients. 


C.2. 3  S  Coefficients. 


(Ti  =  4-8ii  + 


2-8. 


41 


(T2  —  4-8i2  + 


(T3  —  4-8i3  + 


n 

2-842 

n 

2-843 

n 


Si  —  -821  — 


2-8 


31 


S2  —  -822  ~ 


n 

2-8 


32 


S3  =  B. 


n 

2-833 


23 


n 


(C.8) 


(C.9) 
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Appendix  D.  Translational  Osculation  for  the  J2  Perturbed  Circular 

Chief 

This  appendix  derives  the  trajectory  of  the  translational  osculation  described  in  Sec. 
4.2.2.  Beginning  with  the  expressions  for  Xd  and  yd  as 


Xd  =  ^ - 

n 

2x 

yd  =  y - 

n 


(D,l) 


Using  the  closed-form  non-hyperbolic  expressions  given  in  [37],  the  in-plane  relative 
positions  and  velocities  are  given  as 


X 

cos  9 

lsm9 

y 

^sm9 

9 

cos  9 

X 

—9  sin  9 

^  cos  9 

y 

-^cos9 

L  9 

—9  sin  9 

Xq 


yo 


(D.2) 


where  9  =  gnt  and  theta  =  gn,  with  g  dehned  previously  in  Sec.  4.2.2.  Evaluating  Xd 
and  yd  at  an  epoch  time,  and  inserting  the  SS  relative  velocities,  we  can  express  the 
Xdo  and  ydo  directly  proportional  to  xq  and  yo-  This  is  seen  via  Eq.  D.3 


Xdo  =  4x0  +  2  — 
n 


=  Axn  +  2 


-2cnxr 


n 


=  4  (1  -  c)  xo 

/2xo 

l/dO  —  1/0  ~  2  - 

\  n 


=  1/0-2 


ng'^yo 

2c 

n 


=  1  -  -  h/0 

c 


(D.3) 
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After  some  simplification,  substitution  of  the  initial  conditions  of  Eq.  D.3  into 
Eq.  D.2,  and  the  resulting  substitution  into  Eq.  D.l  yields  the  following  closed  form 
expressions  for  Xd  and  yd 


Xd 

cos  9  K,sm9 

^dO 

Vd 

—K~^sm9  cos  9 

VdO 

(D.4) 


where  k  =  This  looks  very  familiar  when  compared  to  a  parameterized  ellipse. 

We  can  now  attempt  to  quantify  the  geometry  of  this  ellipse.  For  now,  assume  a 
dummy  magnitude  C  such  that  Eq.  D.4  is  written  as 


fxdocosd  KydQiiind\ 

\  C  ^  C  ) 

/  a;dosin6'  ydocos9\ 

V  kC  ^  C  ) 


(D.5) 


Now,  dehne  the  following  dummy  phase  angles  as 


sin  (5 
cos  (5 


^dO 

~c 

K-ydo 

C 


Substituting  Eq.  D.6  in  Eq.  D.5  now  gives 


(D,6) 


Xd  =  C  (sin  5  cos  6  +  cos  6  sin  6) 
yd  =  —Cm  (sin  <5  sin  6^  —  cos  5  cos  6) 

where  m  =  k~^ .  Applying  angle  addition  trig  identities 


(D.7) 


sin  (a  +  /5)  =  sin  a  cos  f3  +  cos  a  sin  f3 
cos  (a  +  /5)  =  cos  a  cos  f3  —  sin  a  sin  f3 


(D,8) 
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we  can  collapse  the  Xd  and  yd  expressions  to 


Xd  =  C  sin  {5  +  9)  yd  =Cmcos{5  +  9)  (D.9) 

We  now  need  to  find  the  magnitude  C  by  examining  the  Xd  expression.  Looking  at 
the  Xd  component  of  Eq.  D.9  and  comparing  with  the  Xd  component  of  Eq.  D.5,  we 
observe  the  equality 


C  sin  (5  +  6^)  =  Xdo  cos  9  +  nydo  sin  9  (D.IO) 

Using  the  expression  for  the  linear  combination  of  trigonometric  functions 

a  sin  a;  +  6  cos  a;  =  r  sin  (x  +  a)  (D.ll) 

where  a  =  sin“^  (b/r)  and  r  =  +  6^.  Defining  the  oscillation  C  as  T,  we  find  the 

expression  for  the  oscillation  as 


Such  that  the  final  expression  for  the  Xd  and  yd  trajectory  is 

Xd  =  T  sin  {9  +  (5o) 
yd  =  Tmcos(6'  +  (5o) 

where 

^  ^ 

9  =  grit 

do  =  sm  j 

-c 


(D.12) 


(D.13) 


(D.14) 
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This  is  the  equation  for  an  ellipse  centered  at  the  origin  with  a  semi-major  axis 
of  T/k  and  a  semi-minor  axis  of  T.  The  value  for  k  is  a  hxed  constant.  The  inverse 
of  this  value  is  approximately  1.5  (1.499574...).  This  implies  that  for  a  circular  chief 
perturbed  by  J2,  the  motion  of  the  center  of  the  2x1  osculating  ellipse  follows  an 
approximately  1.5  x  1  ellipse  with  magnitude  solely  determined  by  the  initial  radial 
and  in-track  displacements. 
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Appendix  E.  Circularizing  the  Virtual  Chief  Parameters 

The  following  is  a  description  of  the  mathematics  to  enforce  the  circular  chief  assump¬ 
tion  into  the  six  geometric  parameters  dehning  the  initial  conditions  of  the  Virtual 
Chief  model.  This  is  not  a  primary  objective  of  this  study  and  is  an  exercise  desired 
by  the  research  sponsor  .  The  process  begins  with  a  description  of  the  parameters, 
simplihcations  with  the  circular  assumption,  and  follows  with  the  algebraic  manipu¬ 
lation  and  observations  to  correlate  the  parameters  with  the  relative  orbit  elements. 
The  six  Virtual  Chief  parameters  are 

•  Ai\  An  indication  of  the  scale  of  the  periodicity  of  the  deputy  motion 

•  (j)i\  Initial  phase  angle  in-plane;  also  controls  deputy’s  epoch  location  along  with 
the  skewness  and  orientation  of  the  trajectory 

•  A2.  A  scaling  of  the  drift  rate;  approximate  in-plane  bounded  center  is  also 
located  at  (0,  A2,  0) 

•  02 :  Directly  controls  the  drift  rate 

•  Zmax'-  Maximum  amplitude  of  the  deputy’s  out  of  plane  motion 

•  Tq:  Initial  phase  angle  of  the  deputy’s  out  of  plane  motion 
Dehning  the  Virtual  Chief  parameters  as  the  set  S  such  that 

=  [Ai0iA202^maa;'ho]  (E-1) 
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The  set  can  then  be  dehned  by  the  following  six  relations 


A,  =  ^^/D^  +  {3A  +  2Cy 
01  =  ncto  +  atan2(3A  +  2C,  D) 
dl2  =  y^iB  -2Dy  +  {AA  +  2C)^ 
02  =  atan2(4A  +  2C,B  -  2D) 


^m.ax  — 


To  =  atan2(2:o,  — ) 
n^ 


Where  ric  is  the  mean  motion  of  the  chief,  and  the  variables  A,  B,C,  and  D  in  Equation 
E.2  are  intermediate  parameters  dehned  as 


A  =  cos(z/c(,  -  Mco)xo  -  sin(z/co 
B  =  sin(z/c(,  -  Mco)xo  +  cos(z/c(, 

C  =  — sin(z/co  -  M^){xq  -  (z>co 

TIq 

D  =  — cos(z/co  -  Mco){xo  -  {l>co 

Tic 

The  variables  u  and  M  in  Equation  E.3  refer  to  the  true  and  mean  anomalies  of 
the  chief,  respectively.  Immediately,  the  set  A,  B,  C,  and  D  can  be  reduced  enforcing 
the  circular  assumption  by  noting  that  for  a  body  in  a  circular  orbit,  the  true  and 
mean  anomalies  will  be  equivalent;  this  immediately  will  force  the  cosine  terms  to 
unity  and  null  the  sine  terms  in  the  parameters.  The  next  simplihcation  comes  from 
noting  that  the  time  derivative  of  the  true  anomaly  is  equivalent  to  the  mean  motion 
of  the  circular  chief.  Substitution  of  these  observations  allow  Equation  E.3  to  become 


Mco)yo 

Mco)yo 

-  nc)yo)  +  cos(z/co  -  Mc^){yo 

-  ric)yo)  +  sin(z/co  -  McJ(?/o 


(«>co  -  nc)xo) 
(«>co  -  nc)xo) 


(E.3) 
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A  =  xo 
B  =  Vo 


C 

D 


Tic 

Xo 


(E.4) 


Uc 

Having  the  circular  intermediate  parameters  (Eq  E.4),  we  can  now  reduce  the  com¬ 
ponents  of  S  individually. 


E.l  Simplifying  Ai 

Enforcing  the  circular  assumption  into  Ai,  the  parameter  becomes 


Hi  =  -^D^  +  {3A  +  2C)^ 


(E.5) 


An  immediate  observation  is  the  similarity  to  the  semi-major  axis  of  the  relative  orbit 
in  the  Clohessy- Wiltshire  equations  from  the  ROE  set, 


ae  =  2^/(-)2  +  (3a;+^)2 

Ur  Ur 


(E.6) 


By  obvious  inspection,  the  hrst  parameter  is  a  linear  mapping  to  the  relative  orbit 
element  Oe  by  a  simple  scaling  factor.  Evaluating  Equation  E.6  at  epoch,  it  is  apparent 
that 


Hi  —  -Oe 


(E.7) 


E.2  Simplifying  01 

Earlier  it  was  noted  that 


01  =  uHo  +  atan2(3H  -|-  20,  D) 


(E.8) 
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Using  the  relation  for  0i,  and  substituting  in  the  circular  constraints,  and  assumping 
that  to  is  equivalently  0,  01  becomes 

01  =  Ucto  +  atan2(3A  +  2U,  D) 

=  atan2(3A  +  2U,  D)  (E-9) 

/  20  X  , 

=  atan2(3a;  H - ,  — ) 

Uc  ric 

The  ROE  parameter  0  is  dehned  as 

0  =  atan2(i;,  “incX  +  20)  (E.IO) 


The  two  expressions  in  Equations  E.9  and  E.IO  are  very  similar.  In  fact,  the 
argument  in  the  atan2  functions  are  simply  inverted,  and  the  following  identity  can 
be  utilized 


arctan(a; 


TT 

2 


arctan(a;) 


(E.ll) 


This  is  a  well  known  triginometric  identity  and  is  used  without  proof,  but  observing  a 
simple  right  triangle  will  yield  this  same  result.  Using  this  expression  and  observing 
the  epoch  conditions  on  0 


01  =  atan2(3a;  H — — ) 

Tie  Uc 

0  =  atan2(  — ,  3a;  H — -) 

fir  fir 


(E.12) 


The  claim  is  made  that  the  parameter  0i  condenses  to  the  following  for  a  circular 
chief 


01  = 


(E.13) 
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E.3  Simplifying  A2 


The  parameter  A2  is  defined  in  E.2  and  upon  the  circular  substitution,  the 
following  is  found 

A2  =  -  2D)2  +  (4A  +  2^)2 

(E.14) 


... 


=  J(p--y  +  (4x  +  2^f 


Hr 


Tlr 


Two  components  of  the  ROE  set  describe  the  instantaneous  center  {xd,yd)  of 
the  relative  ellipse  by  the  following  relation 


/I  ^  ‘^y 

Xd  =  4x  ^ - 

Tie 

2x 


(E.15) 


Vd  =  V 


Ur 


Evaluating  Equation  E.15  at  epoch,  it  is  immediately  observed  that  Xd  and  pd 
are  the  parenthetical  expressions  in  Equation  E.14,  and  A2  conveniently  collapses  to 


^2  -  ^l^d  +  Vd 


(E.16) 


This  implies  that  for  the  circular  chief,  the  Virtual  Chief  parameter  A2  is  a 
direct  measurement  of  the  2-norm  of  the  in-plane  position  vector  from  the  chief  to 
the  center  of  the  relative  ellipse. 


E.4  Simplifying  02 

The  expression  for  02  uses  the  same  analytical  functions  as  the  expression  for 
A2.  From  Equation  E.2, 


02  =  atan2(4R  +  2C,B-  2D) 


(E,17) 


112 


Having  already  found  that  the  expressions  for  4A-h2C  and  B  —  2D  are  equivalent 


to  Xd  and  yd,  the  direct  substitution  is  made  such  that 

02  =  atcm2{xd,yd) 


(E.18) 


This  angle  can  be  interpreted  as  the  angle  from  the  center  of  the  relative  ellipse 
to  the  chief  at  epoch.  The  relation  can  be  observed  in  Figure  E.l 


x(Radial) 


Figure  E.l:  02  Visualization 


E.5  Simplifying  Zmax 


A  direct  observation  between  the  ROE  parameter  z^ax  and  the  Virtual  Chief 
yields  duality.  That  is 


^max  —  ^max 


(E.19) 
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E.6  Simplifying  ^ 

The  final  Virtual  Chief  parameter  to  map  and  simplify  is  T,  the  expression  of 
which  is  given  by 


T  =  atan2{z,  — )  (E.20) 

He 

Multiplying  the  arguments  Equation  E.20  by  unity  in  the  form  of  {nc/fic),  the 
expression  is  equivalent  to 


T  =  atan2(2:nc,  i)  (E.21) 

The  expression  in  Equation  E.21  is,  by  observation,  equivalent  to  the  out  of 
plane  ROE  parameter,  also  termed 


E.  7  Boundedness 

For  a  bounded  relative  orbit,  the  Virtual  Chief  condition  is 

0  =  2A  +  C 

0  =  R2sin(02) 

This  will  hold  valid  for  two  constraints  in  ROE  space 
A2  =  OV02  e  ^ 

02  =  i:k7r(k  =  0, 1, 2,  ...)VR2  £  3^ 


(E.22) 


(E.23) 


A  value  for  02  that  equates  to  an  integer  value  of  tt  must  imply  that  (from 
Equation  E.18)  Xd  is  a  zero  value,  while  yd  can  take  on  any  nonzero  value  (avoiding 
the  discontinuity).  A  zero  value  for  A2  implies  that  =  —y"^,  which  will  not  old  in 
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real  space  unless  both  values  are  equivalently  zero.  The  implication  in  ROE  space 
then  follows  as 


02  =  ±k7i{k  =  0, 1,  2, ...)  ^  {xd,  Vd)  =  (0,  Vd) 
^2  =  0  ^  {xd,  Vd)  =  (0,  O)V02  ^  0  e  3? 


(E.24) 


E.  8  Summary 

To  summarize,  allowing  the  ROE  set  composed  of 


Vd  /3  0  Zfnax 


(E.25) 


to  be  denoted  as  A,  then  the  mapping  from  A  to  H  for  a  circular  chief  is  the  following 


02  =  atan2{xd,yd) 

T  =  T 

^max  —  ^max 


(E.26) 
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while  the  inverse  transformation  is  given  by 


O-e  —  77^1 


'5  =  1- 


^max  ^max 


^  ^ 


Xd  = 


Vd  = 


A2tan(02) 

A/rTtar?(^ 

A2 

a/1  +  tan2(02y 


=  ^2  sin(02) 
=  A2  cos(02) 


(E.27) 


Worth  noting  is  that  the  (xd,  Vd)  transformation  is  a  polar  representation  of  the 
in-plane  motion  as  evident  by  Fig.  E.l.  The  transformation  exists  in  both  the  forward 
and  inverse  directions  without  encountering  singularities. 
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Appendix  F.  Stationary  Orbit  Initialization 

A  popular  idea  in  the  field  of  relative  satellite  formation  is  that  of  formation  main¬ 
tenance.  In  this  respect,  the  term  stationary  orbit  applies  to  a  relative  orbit  whose 
trajectory  repeats  itself  periodically.  Derived  equations  of  motion  can  be  investigated 
with  respect  to  projected  secular  growth  and  constraints  on  initial  conditions  can  be 
made  to  hinder  this.  For  example,  in  the  HCW  derivation  the  condition  ijQ  =  —2nxo 
is  found  by  analyzing  the  y  equation  to  remove  secular  terms.  For  unperturbed  dy¬ 
namics,  this  repeated  orbit  results  from  an  equal  periodicity  condition,  which  in  turn 
results  in  equal  semi-major  axes  between  the  chief  and  deputy.  Maintaining  equal 
orbital  periods  between  the  chief  and  deputy  in  the  perturbed  case  is  a  more  dif¬ 
ficult  task  in  that  the  apsides  rotation  modifies  the  nominal  periods  of  the  orbits. 
This  section  presents  cases  to  bound  the  orbit  using  relative  orbit  elements  for  the 
unperturbed  elliptical  chief,  and  the  J2  perturbed  circular  chief. 


F.  0.0.1  Bounding  the  Unperturbed  Elliptical  Chief.  The  equal  peri¬ 
odicity  constraint  yields  an  equal  semi-major  axis  constraint  on  the  motion  for  a 
bounded  orbit  [Sa  =  0).  Using  the  linearized  mapping  from  the  Hill  frame  to  orbital 
element  differences  provided  in  [30]  and  [8],  the  expression  for  Sa  is  given  as 


6a 


2a  (2  -|-  3^1  -|-  2^2)  X  -|-  2av  (1 


^  ,  2a^vp  . 

2ki  +  K2)yS - —X  + 

Vt 


2a 

Vt 


(1  -|-  2ki  -|-  K2)  y 

(F.l) 
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The  coefficients  a,  Ki,  ^2,  andn  can  be  found  in  [30].  However,  making  the  assumption 
that  the  epoch  condition  is  perigee,  the  coefficients  and  other  parameters  simplify  to 


K-  =  0 

a  fp 

Ki  =  -  I - 1 

r  Vr 

K2=0 


Vt  =  TpU^ 


(1  +  e) 

(1-ef 


a 

a  =  — 

Tp 

rp  =  a{l-e) 


which  reduces  the  Sa  expression  to 


(F.2) 


(5a  =  0 


a  /  alp  \\  2a  f  a  f  p  \\  (F-3) 

=  2-  2  +  3-  -  1  xo  +  -  1  +  2-  -  1  Uo 


Tp  \rp 


V, 


Tp  \rp 


Further  simplihcation  results  in 


If  3e  \  2a  f  2a  f  p 

0  =  -  2  -  - - Uo  +  -  1  +  - r  . 

1  —  e  \  1  ~  6/  Ft  V  a  (1  —  e)  \a  (1  —  e) 


-  1 


(F.4) 


And  this  result  further  simplihes  to 


^  ,  yo  /(1-e)^ 


(F.5) 


Finally,  Eq.  F.5  can  be  written  in  general  as 


Vo 


n  (2  +  e) 


\/(l  +  e)(l-y 


(F.6) 
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which  is  near  exact  to  the  constraint  given  in  [22],  To  convert  this  expression  to  ROE 
space,  we  employ  the  Cartesian  to  ROE  relation  as 


®e0  n  I 

Xo  =  — —  COS  15^  +  Xrfo 


?/0  =  tteoncosPo  - 


3nx. 


do 


(F.7) 


Substituting  the  ROE  expressions  in  Eq.  E.6  and  allowing  (  = 


(2+e) 


V  (l+e)(l-e)3 


tteO  COS  Po  -  ^ 

-'^COSpo  +  XdO 


=  c 


(F.8) 


The  following  derives  the  relationship  between  aeo,Xdo,  and  Pq 


SnXdO  y.  /  CleO  rt  I  ^ 

aeon  cos  To - ^ — a  =  C  ^  cos  To  +  j 

n  C®eO  >  3nXdO 

aeon  cos  To  +  =  (xdo  H - ^ — 

/  CcosTo\  ,  3n 

Oeo  1  n  cos  To  H - ^ —  I  =  a^do  I  C  +  Y 


(F.9) 


At  perigee  (To  =  0)  this  becomes  The  final  constraint  for  bounded  motion  becomes 


fleo  _  f  C  +  \ 

Xdo  V2C  +  3ny 


(F.IO) 


Note  that  for  e  =  0,  this  expression  reduces  to  Xdo  =  0  which  is  known  from  [39]  to 
produce  a  bounded  orbit. 

Figure  F.l  shows  the  behavior  of  this  boundedness  requirement  (here  termed 
the  radial  ratio)  and  (.  In  this  example  a  chief  of  semi-major  axis  8000  km  is  used  to 
evaluate  the  value  of  n.  Instantly,  it  is  observed  that  the  radial  ratio  remains  near  the 
same  constant  value,  and  that  as  the  eccentricity  approaches  unity,  the  value  for  the 
radial  ratio  approaches  0.5.  This  is  easily  observed  by  examining  that  as  eccentricity 
approaches  unity,  the  C,  parameter  approaches  infinity,  and  the  radial  ratio  from  Eq. 
F.IO  approaches  0.5. 
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Chief  Eccentricity 

Figure  F.l:  Boundedness  Parameters  as  Functions  of  Chief  Eccentricity  with  a  Chief 
Semi-major  Axis  of  8000  km 

Also,  Fig.  F.2  provides  a  display  of  the  radial  ratio  as  the  chief  semi-major  axis 
is  varied.  The  general  trend  is  followed  for  each  different  value  of  the  semi-major  axis, 
but  the  initial  conditions  require  an  increase  in  radial  displacement. 


F.0.0.2  Bounding  the  J2  Perturbed  Circular  Chief.  In  bounding  the 
relative  trajectory  of  the  circular  chief  perturbed  by  the  J2  effect  there  are  two  possible 
routes.  In  the  Schweighart-Sedwick  model  [7],  initial  conditions  are  derived  to  avoid 
secular  drift  in  a  purely  linear  sense.  Another  approach  is  applying  Schaub’s  method 
for  initializing  J2  invariant  orbits  [30].  This  study  will  determine  expressions  using 
the  Schweighart  initial  conditions. 

The  initial  conditions  given  in  Schweighart  [37]  to  remove  drift  are 


ng 

2c 


?/o  =  -2cnxo 


(F.ll) 
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Figure  F.2:  Radial  Ratio  at  Perigee  as  a  Function  of  Chief  Eccentricity  with  Varying 
Semi-major  Axis 


where  the  coefficients  have  been  previously  dehned.  The  substitution  for  the  ROEs 
can  be  made  into  Eq.  F.ll.  Focusing  on  the  xo  expression 

^nsin/3o  =  ^  (oeo  sin/3o  +  2/do)  (F-12) 

Dividing  Eq.  F.12  by  n,  multiplying  by  2c,  and  collecting  Oeosin/do  terms,  the  condi¬ 
tion  becomes 

aeosin/?o  =  — ^ (F-13) 

c  +  r 

Now,  focusing  on  the  ?/o  expression 


Oeoncos^o  -  =  -2cn  cos^o  +  (F.14) 

Multiplying  the  right  hand  side  of  Eq 
condition  simplihes  to 

tteO  COS  Po 


F.14  through,  and  collecting  like  terms,  the 


3 -4c 
2  -  2c 


XdO 


(F.15) 
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Squaring  and  summing  the  results  of  Eq.  F.13  and  Eq.  F.15 


(aeoCos/3o)V  (aeosin/3o)^  =  Vdo  (F-16) 

Taking  the  square  root  of  Eq.  F.16, 


^J (aeocos/3o)^  +  (aeosin/3o)^  =  y  vlo  (F-17) 

which  directly  implies  the  necessary  condition  to  bound  the  relative  trajectory 

aeo  =  \J  +  Csi/Jo  (F-18) 


where 


C'i  = 
^2  = 


3-4c 
2 -2c, 


c  +  r 


(F.19) 


The  relation  given  in  the  bounded  expression  shows  a  parabolic  surface  dependent  on 
the  values  for  Xdo  and  Udo-  A  very  interesting  note  is  when  these  initial  conditions  are 
used,  the  osculational  translation  of  the  relative  trajectory  follows  a  near  1.5  x  1  ellipse 
with  oscillations  on  very  low  order  magnitude.  The  derivation  of  this  is  provided  in 
Appendix  D. 
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